--- /dev/null
+/****************************************************************************
+ * 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 <iostream.h>
+ #include <fstream.h>
+
+ #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; i<nentr; i++) {
+ AliCascadeVertex *iovertex=new AliCascadeVertex;
+ branch->SetAddress(&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<<ngood<<'\r';
+ if (ngood==100) {
+ cerr<<"Too many good cascades !\n";
+ break;
+ }
+ }
+ if (!in.eof()) cerr<<"Read error (good_cascades) !\n";
+ } else {
+ /*** generate a file with the "good" cascades ***/
+ cerr<<"Marking good cascades (this will take a while)...\n";
+ ngood=good_cascades(gc,100);
+ ofstream out("good_cascades");
+ if (out) {
+ for (Int_t ngd=0; ngd<ngood; ngd++)
+ out<<gc[ngd].nlab<<' '<<gc[ngd].plab<<' '<<
+ gc[ngd].blab<<' '<<gc[ngd].code<<' '<<
+ gc[ngd].px<<' '<<gc[ngd].py<<' '<<gc[ngd].pz<<' '<<
+ gc[ngd].x <<' '<<gc[ngd].y <<' '<<gc[ngd].z <<endl;
+ } else cerr<<"Can not open file (good_cascades) !\n";
+ out.close();
+ }
+
+ /*** create some histograms ***/
+ TH1F *hp=new TH1F("hp","Angular Resolution",30,-30.,30.); //phi resolution
+ hp->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; i<nentr; i++) {
+ AliCascadeVertex *cascade=(AliCascadeVertex*)carray.UncheckedAt(i);
+ cascade->GetV0labels(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; j<ngood; j++) {
+ if (gc[j].code != cascade->GetPdgCode()) 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: ("<<nlab<<","<<plab<<","<<blab<<")\n";
+ continue;
+ }
+
+ pxg=gc[j].px; pyg=gc[j].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+ Double_t phig=TMath::ATan2(pyg,pxg), phi=TMath::ATan2(py,px);
+ Double_t lamg=TMath::ATan2(gc[j].pz,ptg), lam=TMath::ATan2(pz,pt);
+ hp->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; i<ngood; i++) {
+ if (gc[i].code != code) continue;
+ pxg=gc[i].px; pyg=gc[i].py; ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+ hgood->Fill(ptg);
+ nlab=gc[i].nlab; plab=gc[i].plab; blab=gc[i].blab;
+ if (nlab < 0) continue;
+ cerr<<"Cascade ("<<nlab<<','<<plab<<","<<blab<<") has not been found !\n";
+ }
+
+ carray.Delete();
+
+ Stat_t ng=hgood->GetEntries();
+ Stat_t nf=hfound->GetEntries();
+
+ cerr<<"Number of found cascades: "<<nentr<<" ("<<nf<<" in the peak)\n";
+ cerr<<"Number of good cascades: "<<ng<<endl;
+
+ if (ng!=0)
+ cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
+
+ gStyle->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<<np<<'\r';
+ TParticle *cp=gAlice->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; i<ngt; i++) if (label[i]==blab) break;
+ if (i==ngt) continue;
+
+ /** is the V0 "good" ? **/
+ Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
+ if (nlab==plab) continue;
+ if (nlab<0) continue;
+ if (plab<0) continue;
+
+ for (i=0; i<ngt; i++) if (label[i]==nlab) break;
+ if (i==ngt) continue;
+ for (i=0; i<ngt; i++) if (label[i]==plab) break;
+ if (i==ngt) continue;
+
+ /*** fiducial volume ***/
+ Double_t x=bp->Vx(), y=bp->Vy(), z=bp->Vz(), r2=x*x+y*y; //bachelor
+ if (r2<r2min) continue;
+ if (r2>r2max) continue;
+ TParticle *pp=gAlice->Particle(plab);
+ x=pp->Vx(); y=pp->Vy(); r2=x*x+y*y; //V0
+ if (r2<r2min) continue;
+ if (r2>r2max) 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;
+}
--- /dev/null
+#ifndef __CINT__
+ #include <iostream.h>
+ #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;
+}
--- /dev/null
+/**************************************************************************
+ * 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 <iostream.h>
+#include <TMath.h>
+
+#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;
+}
--- /dev/null
+#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 <TObject.h>
+#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
+
+
--- /dev/null
+/**************************************************************************
+ * 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 <TFile.h>
+#include <TTree.h>
+#include <TObjArray.h>
+#include <iostream.h>
+
+#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; i<nentrV0; i++) {
+
+ AliV0vertex *ioVertex=new AliV0vertex;
+ branch->SetAddress(&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; i<nentr; i++) {
+
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ branch->SetAddress(&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; i<nV0; i++) {
+
+ AliV0vertex *V0ver=(AliV0vertex *)vtxV0.UncheckedAt(i);
+
+ V0ver->ChangeMassHypothesis(kLambda0); //I.B.
+
+ if (V0ver->GetEffMass()<massLambda-fMassWin || // condition of the V0 mass window (cut fMassWin)
+ V0ver->GetEffMass()>massLambda+fMassWin) continue;
+
+ if (V0ver->GetD(0,0,0)<fDV0min) continue; // condition of minimum impact parameter of the V0 (cut fDV0min)
+ // here why not cuting on pointing angle ???
+
+ // for each vertex in the good mass range, loop on all tracks (= bachelor candidates)
+
+ for (Int_t k=0; k<ntrack; k++) {
+
+ AliITStrackV2 *bachtrk=(AliITStrackV2 *)trks.UncheckedAt(k);
+
+ if (TMath::Abs(bachtrk->GetD())<fDBachMin) continue; // eliminate to small impact parameters
+
+ if (V0ver->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 (cost<fCPAmax) continue; // condition on the cascade pointing angle
+
+ //cascade.ChangeMassHypothesis(); //default is Xi
+
+
+ ioCascade=&cascade; cascTree.Fill();
+
+ ncasc++;
+
+ }
+ }
+
+ cerr<<"Number of reconstructed cascades: "<<ncasc<<endl;
+
+ cascTree.Write();
+
+ trks.Delete();
+ vtxV0.Delete();
+
+ savedir->cd();
+
+ 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;
+}
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+#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
+
+
cerr<<"Track "<<lab<<" was not found !\n";
continue;
}
- 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);
#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 = "<<SlowOrFast<<endl;
-///////////////// Conversion AliITSRecPoint -> AliITSclusterV2 //////////////
-
- // for the fast recpoints
- if(SlowOrFast=='f') {
- cout<<"22 !!!! SlowOrFast = "<<SlowOrFast<<endl;
- /*
- if (gAlice) {delete gAlice; gAlice=0;}
-
- TFile *in=TFile::Open("galice.root","update");
- 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 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);
cerr<<"Number of entries: "<<nentr<<endl;
- Float_t kmip; // ADC->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; i<nentr; i++) {
points->Clear();
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 ="<<i<<","<<ncl<<endl;
+ Float_t kmip=1; // ADC->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; j<ncl; j++) {
AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(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="<<lab[0]<<" "<<endl;
-
- if(label>=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: "<<part->P()<<endl; break;}
- label=m;
- part=(TParticle*)gAlice->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: "<<part->P()<<endl; break;}
+ label=m;
+ part=(TParticle*)gAlice->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);
}
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;
-Int_t AliITStestV2() {
+Int_t AliITStestV2(Char_t SlowOrFast='s') {
Int_t rc=0;
if (gAlice) {delete gAlice; gAlice=0;}
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;
// 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 <TMatrixD.h>
fC40=t.fC40; fC41=t.fC41; fC42=t.fC42; fC43=t.fC43; fC44=t.fC44;
Int_t n=GetNumberOfClusters();
- //for (Int_t i=0; i<n; i++) fIndex[i]=t.fIndex[i];
- //b.b.
for (Int_t i=0; i<n; i++) {
- fIndex[i]=t.fIndex[i];
- fdEdxSample[i]=t.fdEdxSample[i];
+ fIndex[i]=t.fIndex[i];
+ if (i<4) fdEdxSample[i]=t.fdEdxSample[i];
}
-
}
-/*
+
//_____________________________________________________________________________
Int_t AliITStrackV2::Compare(const TObject *o) const {
//-----------------------------------------------------------------
// This function compares tracks according to the their curvature
//-----------------------------------------------------------------
AliITStrackV2 *t=(AliITStrackV2*)o;
- Double_t co=TMath::Abs(t->Get1Pt());
- 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 (c<co) return -1;
return 0;
}
-*/
-//_____________________________________________________________________________
-Int_t AliITStrackV2::Compare(const TObject *o) const {
- //-----------------------------------------------------------------
- // This function compares tracks according to the their curvature
- //-----------------------------------------------------------------
- AliITStrackV2 *t=(AliITStrackV2*)o;
-
- Double_t p2=1./(Get1Pt()*Get1Pt());
- Double_t b2=p2/(p2 + GetMass()*GetMass());
- Double_t po2=1./(t->Get1Pt()*t->Get1Pt());
- Double_t bo2=po2/(po2 + t->GetMass()*t->GetMass());
- if (p2*b2>po2*bo2) return -1;
- else if (p2*b2<po2*bo2) return 1;
- return 0;
-}
//_____________________________________________________________________________
void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
}
//____________________________________________________________________________
-Int_t AliITStrackV2::PropagateToVertex(Double_t x0,Double_t rho) {
+Int_t AliITStrackV2::PropagateToVertex(Double_t d,Double_t x0) {
//------------------------------------------------------------------
//This function propagates a track to the minimal distance from the origin
//------------------------------------------------------------------
Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
- Propagate(fAlpha,xv,0.,0.);
+ PropagateTo(xv,d,x0);
return 0;
}
//------------------------------------------------------------------
Double_t dx=xk-fX;
Double_t f1=fP2, f2=f1 + fP4*dx;
- if (TMath::Abs(f2) >= 0.99999) {
+ if (TMath::Abs(f2) >= 0.9999) {
Int_t n=GetNumberOfClusters();
if (n>kWARN)
cerr<<n<<" AliITStrackV2::GetGlobalXYZat: Propagation failed !\n";
}
//____________________________________________________________________________
-Int_t
-AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho) {
+Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) {
+ //------------------------------------------------------------------
+ //This function corrects the track parameters for crossed material
+ //------------------------------------------------------------------
+ Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+ Double_t beta2=p2/(p2 + GetMass()*GetMass());
+ d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+
+ //Multiple scattering******************
+ if (d!=0) {
+ //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
+ Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*d*9.36*2.33;
+ fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
+ fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
+ fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
+ fC44 += theta2*fP3*fP4*fP3*fP4;
+ }
+
+ //Energy losses************************
+ if (x0!=0.) {
+ d*=x0;
+ Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
+ fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
+ }
+
+ if (!Invariant()) return 0;
+
+ return 1;
+}
+
+//____________________________________________________________________________
+Int_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
//------------------------------------------------------------------
//This function propagates a track
//------------------------------------------------------------------
Double_t x1=fX, x2=xk, dx=x2-x1;
Double_t f1=fP2, f2=f1 + fP4*dx;
- if (TMath::Abs(f2) >= 0.99999) {
+ if (TMath::Abs(f2) >= 0.9999) {
Int_t n=GetNumberOfClusters();
if (n>kWARN)
cerr<<n<<" AliITStrackV2::PropagateTo: Propagation failed !\n";
fX=x2;
- Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + GetMass()*GetMass());
-
- //Multiple scattering******************
- //x0=0.;
- if (x0!=0) {
- x0*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
- Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
- fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
- fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
- fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
- fC44 += theta2*fP3*fP4*fP3*fP4;
- }
-
- //Energy losses************************
- if (rho!=0.) {
- rho*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
- Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*rho;
- if (x1 < x2) dE=-dE;
- fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
- }
-
- if (!Invariant()) {cout<<"Propagate !\n"; return 0;}
+ if (!CorrectForMaterial(d,x0)) return 0;
return 1;
}
return 1;
}
-
-//____________________________________________________________________________
-Int_t AliITStrackV2::Update(const Double_t* m, Double_t chi2, UInt_t index) {
- //------------------------------------------------------------------
- //This function updates track parameters with a vertex constraint
- //------------------------------------------------------------------
- Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
- Double_t c00=fC00;
- Double_t c10=fC10, c11=fC11;
- Double_t c20=fC20, c21=fC21, c22=fC22;
- Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
- Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
-
-
- TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
- TMatrixD C(5,5);
- C(0,0)=fC00;
- C(1,0)=fC10; C(1,1)=fC11;
- C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
- C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
- C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
-
- C(0,1)=C(1,0);
- C(0,2)=C(2,0); C(1,2)=C(2,1);
- C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
- C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
-
- TMatrixD H(4,5); H.UnitMatrix();
- TMatrixD Ht(TMatrixD::kTransposed,H);
- TVectorD mm(4); mm(0)=m[0]; mm(1)=m[1]; mm(2)=m[2]; mm(3)=m[3];
- TMatrixD V(4,4);
- V(0,0)=m[4 ]; V(0,1)=m[5 ]; V(0,2)=m[6 ]; V(0,3)=m[7 ];
- V(1,0)=m[8 ]; V(1,1)=m[9 ]; V(1,2)=m[10]; V(1,3)=m[11];
- V(2,0)=m[12]; V(2,1)=m[13]; V(2,2)=m[14]; V(2,3)=m[15];
- V(3,0)=m[16]; V(3,1)=m[17]; V(3,2)=m[18]; V(3,3)=m[19];
-
- TMatrixD tmp(H,TMatrixD::kMult,C);
- TMatrixD R(tmp,TMatrixD::kMult,Ht); R+=V;
-
- R.Invert();
-
- TMatrixD K(C,TMatrixD::kMult,Ht); K*=R;
-
- TVectorD savex=x;
- x*=H; x-=mm; x*=-1; x*=K; x+=savex;
-
- TMatrixD saveC=C;
- C.Mult(K,tmp); C-=saveC; C*=-1;
-
- fP0=x(0); fP1=x(1); fP2=x(2); fP3=x(3); fP4=x(4);
- fC00=C(0,0);
- fC10=C(1,0); fC11=C(1,1);
- fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
- fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
- fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
-
-
- if (!Invariant()) {
- fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
- fC00=c00;
- fC10=c10; fC11=c11;
- fC20=c20; fC21=c21; fC22=c22;
- fC30=c30; fC31=c31; fC32=c32; fC33=c33;
- fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
- return 0;
- }
-
- Int_t n=GetNumberOfClusters();
- fIndex[n]=index;
- SetNumberOfClusters(n+1);
- SetChi2(GetChi2()+chi2);
-
- return 1;
-}
-
Int_t AliITStrackV2::Invariant() const {
//------------------------------------------------------------------
// This function is for debugging purpose only
//------------------------------------------------------------------
Int_t n=GetNumberOfClusters();
- //if (TMath::Abs(fP1)>11.5)
- //if (fP1*fP4<0) {
- // if (n>kWARN) cerr<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
-
- if (TMath::Abs(fP2)>=1) {if (n>kWARN) cerr<<"fP2="<<fP2<<endl; return 0;}
-
- if (fC00<=0) {if (n>kWARN) cerr<<"fC00="<<fC00<<endl; return 0;}
- if (fC11<=0) {if (n>kWARN) cerr<<"fC11="<<fC11<<endl; return 0;}
- if (fC22<=0) {if (n>kWARN) cerr<<"fC22="<<fC22<<endl; return 0;}
- if (fC33<=0) {if (n>kWARN) cerr<<"fC33="<<fC33<<endl; return 0;}
- if (fC44<=0) {if (n>kWARN) cerr<<"fC44="<<fC44<<endl; return 0;}
- /*
- TMatrixD m(5,5);
- m(0,0)=fC00;
- m(1,0)=fC10; m(1,1)=fC11;
- m(2,0)=fC20; m(2,1)=fC21; m(2,2)=fC22;
- m(3,0)=fC30; m(3,1)=fC31; m(3,2)=fC32; m(3,3)=fC33;
- m(4,0)=fC40; m(4,1)=fC41; m(4,2)=fC42; m(4,3)=fC43; m(4,4)=fC44;
-
- m(0,1)=m(1,0);
- m(0,2)=m(2,0); m(1,2)=m(2,1);
- m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
- m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
-
- Double_t det=m.Determinant();
-
- if (det <= 0) {
- if (n>kWARN) { cerr<<" bad determinant "<<det<<endl; m.Print(); }
- return 0;
+ if (TMath::Abs(fP2)>=0.9999){
+ if (n>kWARN) cout<<"AliITStrackV2::Invariant : fP2="<<fP2<<endl;
+ return 0;
+ }
+ if (fC00<=0 || fC00>9.) {
+ if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC00="<<fC00<<endl;
+ return 0;
+ }
+ if (fC11<=0 || fC11>9.) {
+ if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC11="<<fC11<<endl;
+ return 0;
+ }
+ if (fC22<=0 || fC22>1.) {
+ if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC22="<<fC22<<endl;
+ return 0;
+ }
+ if (fC33<=0 || fC33>1.) {
+ if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC33="<<fC33<<endl;
+ return 0;
+ }
+ if (fC44<=0 || fC44>6e-5) {
+ if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC44="<<fC44<<endl;
+ return 0;
}
- */
return 1;
}
//____________________________________________________________________________
-Int_t
-AliITStrackV2::Propagate(Double_t alp,Double_t xk,Double_t x0,Double_t rho) {
+Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
//------------------------------------------------------------------
//This function propagates a track
//------------------------------------------------------------------
Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
Double_t pp2=fP2*ca - cf*sa;
- if (TMath::Abs(pp2) >= 0.99999) {
+ if (TMath::Abs(pp2) >= 0.9999) {
Int_t n=GetNumberOfClusters();
if (n>kWARN)
cerr<<n<<" AliITStrackV2::Propagate: Rotation failed !\n";
cf=ca + sf*sa/cf;
- if (!Invariant()) {cout<<dalp<<" Rotate !\n"; return 0;}
+ if (!Invariant()) return 0;
x1=fX; Double_t x2=xk, dx=x2-x1;
Double_t f1=fP2, f2=f1 + fP4*dx;
- if (TMath::Abs(f2) >= 0.99999) {
+ if (TMath::Abs(f2) >= 0.9999) {
Int_t n=GetNumberOfClusters();
if (n>kWARN)
cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
- pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
-
- //Multiple scattering******************
- //x0=0.;
- if (x0!=0.) {
- x0*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
- Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
- fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
- fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
- fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
- fC44 += theta2*fP3*fP4*fP3*fP4;
- }
-
- //Energy losses************************
- if (rho!=0.) {
- rho*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
- Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*rho;
- if (x1 < x2) dE=-dE;
- fP4*=(1.- sqrt(pp2+GetMass()*GetMass())/pp2*dE);
- }
-
if (!Invariant()) {
fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
fC00=c00;
Double_t a=2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2))-fP4*(fX*fX + fP0*fP0);
if (fP4<0) a=-a;
return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
-}
+}
+Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
+ //------------------------------------------------------------------
+ //This function improves angular track parameters
+ //------------------------------------------------------------------
+ Double_t dy=fP0-yv, dz=fP1-zv;
+ Double_t r2=fX*fX+dy*dy;
+ Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
+ Double_t beta2=p2/(p2 + GetMass()*GetMass());
+ x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
+ //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
+ Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
+ {
+ Double_t parp=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
+ Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
+ sigma2p += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
+ sigma2p += kSigmaYV*kSigmaYV/r2;
+ sigma2p += 0.25*fC44*fX*fX;
+ Double_t eps2p=sigma2p/(fC22+sigma2p);
+ fP0 += fC20/(fC22+sigma2p)*(parp-fP2);
+ fP2 = eps2p*fP2 + (1-eps2p)*parp;
+ fC22 *= eps2p;
+ fC20 *= eps2p;
+ }
+ {
+ Double_t parl=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
+ Double_t sigma2l=theta2;
+ sigma2l += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
+ sigma2l += kSigmaZV*kSigmaZV/r2;
+ Double_t eps2l=sigma2l/(fC33+sigma2l);
+ fP1 += fC31/(fC33+sigma2l)*(parl-fP3);
+ fP4 += fC43/(fC33+sigma2l)*(parl-fP3);
+ fP3 = eps2l*fP3 + (1-eps2l)*parl;
+ fC33 *= eps2l; fC43 *= eps2l;
+ fC31 *= eps2l;
+ }
+ if (!Invariant()) return 0;
+ return 1;
+}
+/*
Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
//------------------------------------------------------------------
//This function improves angular track parameters
Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
Double_t beta2=p2/(p2 + GetMass()*GetMass());
x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
- Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
+ //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
+ Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
if (!Invariant()) return 0;
return 1;
}
-
-/*
-Int_t AliITStrackV2::Improve(Double_t x0,Double_t xv,Double_t yv) {
- //------------------------------------------------------------------
- //This function improves angular track parameters
- //------------------------------------------------------------------
- TMatrixD I(5,5);
- TVectorD v(5); v(0)=fP0; v(1)=fP1; v(2)=fP2; v(3)=fP3; v(4)=fP4;
-
- Double_t r2=fX*fX+fP0*fP0;
- Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + GetMass()*GetMass());
- x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
- Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
-
- v(2)=0.5*(fP4*fX + fP0*TMath::Sqrt(4/r2-fP4*fP4));
- Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
- sigma2 += fC00/r2*(1.- fP0*fP0/r2)*(1.- fP0*fP0/r2);
- sigma2 += kSigmaYV*kSigmaYV/r2;
- I(2,2)=1/sigma2;
-
- v(3)=0.5*fP4*fP1/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
- sigma2=theta2;
- sigma2 += fC11/r2+fC00*fP0*fP0*fP1*fP1/(r2*r2*r2);
- sigma2 += kSigmaZV*kSigmaZV/r2;
- I(3,3)=1/sigma2;
-
- Double_t tgl=fP3;
-
- TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
- TMatrixD C(5,5);
- C(0,0)=fC00;
- C(1,0)=fC10; C(1,1)=fC11;
- C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
- C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
- C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
-
- C(0,1)=C(1,0);
- C(0,2)=C(2,0); C(1,2)=C(2,1);
- C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
- C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
-
- TMatrixD tmp(I,TMatrixD::kMult,C),U(5,5); U.UnitMatrix();
- U+=tmp;
- U.Invert();
- TMatrixD W1(U);
- TMatrixD W2(tmp,TMatrixD::kMult,W1);
-
- v*=W2; x*=W1; x+=v;
-
- C*=W1;
-
-
- fP0=x(0); fP1=x(1); fP2=x(2); fP3=x(3); fP4=x(4);
- fC00=C(0,0);
- fC10=C(1,0); fC11=C(1,1);
- fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
- fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
- fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
-
- eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
- fP4*=eps;
- fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
-
- if (!Invariant()) return 0;
- return 1;
-}
*/
-
void AliITStrackV2::ResetCovariance() {
//------------------------------------------------------------------
//This function makes a track forget its history :)
void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
//-----------------------------------------------------------------
- // This funtion calculates dE/dX within the "low" and "up" cuts.
+ // This function calculates dE/dX within the "low" and "up" cuts.
+ // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
//-----------------------------------------------------------------
Int_t i;
- Int_t nc=GetNumberOfClusters();
- if(nc != 6) cout<<"!!!Warning: ncl isn't 6, ="<<nc<<endl;
-
+ Int_t nc=4;
// The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
// Take only SSD and SDD
- nc=4;
- Int_t swap;//stupid sorting
-
- cout<<"Start sorting (low,up,nc)..."<<low<<" "<<up<<" "<<nc<<endl;
-
- /*
- for (i=0; i<6; i++) { // b.b.
- cout<<"! cl befor sort: cl,dEdx ="<<i<<","<<fdEdxSample[i]<<endl;
- }
- */
+ Int_t swap;//stupid sorting
do {
swap=0;
for (i=0; i<nc-1; i++) {
}
} while (swap);
- for (i=0; i<nc; i++) { // b.b.
- cout<<" i, sorted dEdx ="<<i<<","<<fdEdxSample[i]<<endl;
- }
- Int_t nl=Int_t(low*nc), nu=Int_t(up*(nc-1)); //b.b. to take two lowest dEdX
- // values from four ones choose
- // nu=2
-
- cout<<" Cook: nl,nu, dEdX samples ="<<nl<<" "<<nu<<" ";
+ Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
+ // values from four ones choose
+ // nu=2
Float_t dedx=0;
- for (i=nl; i<nu; i++){
- dedx += fdEdxSample[i]; cout<<" "<<fdEdxSample[i]<<" ";}
- cout<<endl;
- dedx /= float(nu);
+ for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
+ dedx /= (nu-nl);
- cout<<"! CookdEdx end: dedx ="<<dedx<<endl;
SetdEdx(dedx);
}
-
-
// ITS Track Class
//
// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
//-------------------------------------------------------------------------
AliITStrackV2():AliKalmanTrack(){}
AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *);
AliITStrackV2(const AliITStrackV2& t);
- Int_t PropagateToVertex(Double_t x0=36.66,Double_t rho=1.2e-3);
- Int_t Propagate(Double_t alpha, Double_t xr, Double_t x0, Double_t rho);
- Int_t PropagateTo(Double_t xr,Double_t x0=21.82,Double_t rho=2.33);
- Int_t Update(const Double_t *m, Double_t chi2, UInt_t i);
+ Int_t PropagateToVertex(Double_t d=1.2e-3*3., Double_t x0=36.66);
+ Int_t Propagate(Double_t alpha, Double_t xr);
+ Int_t CorrectForMaterial(Double_t d, Double_t x0=21.82);
+ Int_t PropagateTo(Double_t xr, Double_t d, Double_t x0=21.82);
Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i);
Int_t Improve(Double_t x0,Double_t yv,Double_t zv);
void SetdEdx(Double_t dedx) {fdEdx=dedx;}
- void SetSampledEdx(Float_t q, Int_t i); // b.b.
- //void CookdEdx(Double_t low=0., Double_t up=1.) {}
- void CookdEdx(Double_t low=0., Double_t up=0.8); // b.b.
+ void SetSampledEdx(Float_t q, Int_t i);
+ void CookdEdx(Double_t low=0., Double_t up=0.51);
void SetDetectorIndex(Int_t i) {SetLabel(i);}
void ResetCovariance();
void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
UInt_t fIndex[kMaxLayer]; // indices of associated clusters
- Float_t fdEdxSample[8]; // array of dE/dx samples b.b.
+ Float_t fdEdxSample[4]; // array of dE/dx samples b.b.
ClassDef(AliITStrackV2,1) //ITS reconstructed track
};
xr=fX;
x[0]=GetY(); x[1]=GetZ(); x[2]=GetSnp(); x[3]=GetTgl(); x[4]=Get1Pt();
}
-//b.b.
+
inline
void AliITStrackV2::SetSampledEdx(Float_t q, Int_t i) {
//----------------------------------------------------------------------
- //
+ // This function stores dEdx sample corrected for the track segment length
+ // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
//----------------------------------------------------------------------
+ if (i<0) return;
+ if (i>3) return;
Double_t s=GetSnp(), t=GetTgl();
- //cout<<"before corr!! i,q="<<i<<" "<<q<<" "<<endl;
q *= TMath::Sqrt((1-s*s)/(1+t*t));
fdEdxSample[i]=q;
- //cout<<"after corr!! i,q="<<i<<" "<<q<<" "<<endl;
}
-
-
-
#endif
// Implementation of the ITS tracker class
//
// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
+// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
//-------------------------------------------------------------------------
#include <TFile.h>
#include <TTree.h>
//#define DEBUG
#ifdef DEBUG
-Int_t LAB=7;
+Int_t LAB=5126;
#endif
AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
nentr=(Int_t)tpcTree->GetEntries();
for (Int_t i=0; i<nentr; i++) {
tpcTree->GetEvent(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;
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<<tpcLabel<<" *****************\n";
#endif
+
try {
ResetTrackToFollow(*t);
} catch (const Char_t *msg) {
}
ResetBestTrack();
- Double_t r2=fTrackToFollow.GetX()*fTrackToFollow.GetX();
- Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2;
- if (constraint) fTrackToFollow.Improve(x0,GetY(),GetZ());
-
- Double_t xk=80.;
- fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there
-
- xk-=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk-=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk-=2.0;
- fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
- xk-=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk-=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
-
- xk=61.;
- fTrackToFollow.PropagateTo(xk,0.,0.); //C02
-
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk -=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk -=0.5;
- fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029); //Nomex
- xk -=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk -=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
-
-
for (FollowProlongation(); fI<kMaxLayer; fI++) {
while (TakeNextProlongation()) FollowProlongation();
}
if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z))
throw "AliITStrackerV2::PropagateBack: failed to estimate track !";
Double_t phi=TMath::ATan2(y,x);
- Double_t d=layer.GetThickness(phi,z);
Int_t idet=layer.FindDetectorIndex(phi,z);
if (idet<0)
throw "AliITStrackerV2::PropagateBack: failed to find a detector !\n";
const AliITSdetector &det=layer.GetDetector(idet);
r=det.GetR(); phi=det.GetPhi();
- if (!fTrackToFollow.Propagate(phi,r,d/21.82*2.33,d*2.33)) throw "";
+ if (!fTrackToFollow.Propagate(phi,r)) throw "";
+ fTrackToFollow.SetDetectorIndex(idet);
const AliITSclusterV2 *cl=0;
Int_t index=0;
Double_t maxchi2=kMaxChi2;
if (l==fI) {
- if (idet != c->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 (chi2<kMaxChi2) {
const AliITSclusterV2 *cc=0; Int_t ci;
while ((cc=layer.GetNextCluster(ci))!=0) {
- if (idet != cc->GetDetectorIndex()) continue;
+ idet=cc->GetDetectorIndex();
+ if (idet != fTrackToFollow.GetDetectorIndex()) continue;
Double_t chi2=fTrackToFollow.GetPredictedChi2(cc);
if (chi2<maxchi2) {
cl=cc; index=(fI<<28)+ci; maxchi2=chi2;
if (cl) {
if (!fTrackToFollow.Update(cl,maxchi2,index))
cerr<<"AliITStrackerV2::PropagateBack: filtering failed !\n";
- continue;
}
+
+ x=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ());
+ fTrackToFollow.CorrectForMaterial(-x);
+
}
- Double_t xk=61.;
- fTrackToFollow.PropagateTo(xk,0.,0.); //Air
-
- xk +=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
- xk +=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk +=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk +=0.5;
- fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029); //Nomex
- xk +=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk +=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk +=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7); //Al
-
- xk=80.;
- fTrackToFollow.PropagateTo(xk,0.,0.); //CO2
-
- xk+=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
- xk+=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk+=2.0;
- fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
- xk+=0.02;
- fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45); //Kevlar
- xk+=0.005;
- fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
+ Double_t xk=52.,x,y,z; fTrackToFollow.GetGlobalXYZat(xk,x,y,z);
+ if (TMath::Abs(y)<7.77)
+ fTrackToFollow.PropagateTo(xk,-0.19,24.);
+ fTrackToFollow.PropagateTo(61,-0.0110);
+ fTrackToFollow.PropagateTo(80.,-0.0053);
fTrackToFollow.SetLabel(itsLabel);
otrack=new AliTPCtrack(fTrackToFollow,fTrackToFollow.GetAlpha());
return 0;
}
+
AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
//--------------------------------------------------------------------
// Return pointer to a given cluster
Int_t tryAgain=kLayersToSkip;
while (fI) {
- fI--;
+ Int_t i=fI-1;
#ifdef DEBUG
-cout<<fI<<' ';
+cout<<i<<' ';
#endif
- AliITSlayer &layer=fLayers[fI];
- AliITStrackV2 &track=fTracks[fI];
+ AliITSlayer &layer=fLayers[i];
+ AliITStrackV2 &track=fTracks[i];
Double_t r=layer.GetR();
- if (fI==3 || fI==1) {
- Double_t rs=0.5*(fLayers[fI+1].GetR() + r);
- Double_t ds=0.034; if (fI==3) ds=0.039;
- Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
- fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,dx0r,dr);
+ if (i==3 || i==1) {
+ Double_t rs=0.5*(fLayers[i+1].GetR() + r);
+ Double_t d=0.011; if (i==3) d=0.0053;
+ if (!fTrackToFollow.PropagateTo(rs,d)) {
+ //cerr<<"AliITStrackerV2::FollowProlongation: "
+ //"propagation failed !\n";
+ break;
+ }
}
//find intersection
if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
//cerr<<"AliITStrackerV2::FollowProlongation: "
//"failed to estimate track !\n";
- break;
+ break;
}
Double_t phi=TMath::ATan2(y,x);
- Double_t d=layer.GetThickness(phi,z);
Int_t idet=layer.FindDetectorIndex(phi,z);
if (idet<0) {
//cerr<<"AliITStrackerV2::FollowProlongation: "
//propagate to the intersection
const AliITSdetector &det=layer.GetDetector(idet);
phi=det.GetPhi();
- if (!fTrackToFollow.Propagate(phi,det.GetR(),d/21.82*2.33,d*2.33)) {
+ if (!fTrackToFollow.Propagate(phi,det.GetR())) {
//cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
break;
}
- if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > 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;
}
}
- if (fI) fI++;
}
-
Int_t AliITStrackerV2::TakeNextProlongation() {
//--------------------------------------------------------------------
- //This function takes another track prolongation
+ // This function takes another track prolongation
+ //
+ // dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
//--------------------------------------------------------------------
- //Double_t m[20];
- Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
-
AliITSlayer &layer=fLayers[fI];
AliITStrackV2 &t=fTracks[fI];
- Int_t &constraint=fConstraint[fPass];
Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
const AliITSclusterV2 *c=0; Int_t ci=-1;
Double_t chi2=12345.;
while ((c=layer.GetNextCluster(ci))!=0) {
-
-#ifdef DEBUG
-//if (fI==0)
-//if (c->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<<fI<<" change detector !\n";
if (TMath::Abs(t.GetZ() - c->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 (chi2<kMaxChi2) break;
+ chi2=t.GetPredictedChi2(c); if (chi2<kMaxChi2) break;
}
#ifdef DEBUG
ResetTrackToFollow(t);
- //if (!fTrackToFollow.Update(m,chi2,(fI<<28)+ci)) {
if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
//cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
return 0;
}
- // b.b. change for the dEdX
- fTrackToFollow.SetSampledEdx(c->GetQ(),
- 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="<<c->GetLabel(0)<<' '
<<fTrackToFollow.GetNumberOfClusters()<<' '
- <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<endl<<endl;
+ <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<' '
+ <<fTrackToFollow.Get1Pt()<<endl<<endl;
#endif
return 1;
}
-
-
AliITStrackerV2::AliITSlayer::AliITSlayer() {
//--------------------------------------------------------------------
//default AliITSlayer constructor
//This function adds a cluster to this layer
//--------------------------------------------------------------------
if (fN==kMaxClusterPerLayer) {
- cerr<<"AliITStrackerV2::AliITSlayer::InsertCluster(): Too many clusters !\n";
- return 1;
+ cerr<<"AliITStrackerV2::AliITSlayer::InsertCluster(): "
+ "Too many clusters !\n";
+ return 1;
}
if (fN==0) {fClusters[fN++]=c; return 0;}
}
Double_t
-AliITStrackerV2::AliITSlayer::GetThickness(Double_t phi, Double_t z) const {
- //--------------------------------------------------------------------
- //This function returns the thickness of this layer
- //--------------------------------------------------------------------
- //-pi<phi<+pi
- if (3 <fR&&fR<8 ) return 2.5*0.096;
- if (13<fR&&fR<26) return 1.1*0.088;
- if (37<fR&&fR<41) return 1.1*0.085;
- return 1.1*0.081;
-}
+AliITStrackerV2::AliITSlayer::GetThickness(Double_t y, Double_t z) const {
+ //--------------------------------------------------------------------
+ //This function returns the layer thickness at this point (units X0)
+ //--------------------------------------------------------------------
+ Double_t d=0.0085;
+
+ if (43<fR&&fR<45) { //SSD2
+ d=0.0036;
+ if (TMath::Abs(y-0.00)>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 (37<fR&&fR<41) { //SSD1
+ d=0.0036;
+ if (TMath::Abs(y-0.00)>3.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 (13<fR&&fR<26) { //SDD
+ d=0.0034;
+ if (TMath::Abs(y-0.00)>3.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 (6<fR&&fR<8) { //SPD2
+ d=0.0093;
+ if (TMath::Abs(y-3.08)>0.45) d+=0.0064;
+ if (TMath::Abs(y-3.03)<0.10) d+=0.0192;
+ } else
+ if (3<fR&&fR<5) { //SPD1
+ d=0.0093;
+ if (TMath::Abs(y+0.21)>0.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)
//--------------------------------------------------------------------
- //-pi<phi<+pi
- Double_t d=0.1;
+ Double_t d=0.1/65.19*1.848;
Double_t xn=fLayers[fI].GetR();
for (Int_t i=0; i<fI; i++) {
Double_t xi=fLayers[i].GetR();
- d+=fLayers[i].GetThickness(phi,z)*xi*xi;
+ d+=fLayers[i].GetThickness(y,z)*xi*xi;
}
if (fI>1) {
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);
}
for (Int_t i=fI; i<fN; i++) {
const AliITSclusterV2 *c=fClusters[i];
if (c->GetZ() > fZmax) break;
- //if (c->IsUsed()) continue;
+ if (c->IsUsed()) continue;
const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
Double_t y=fR*det.GetPhi() + c->GetY();
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:
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() {
+/****************************************************************************
+ * 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 <iostream.h>
#include <fstream.h>
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.;
Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y;
if (r2<r2min) continue;
if (r2>r2max) 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();
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
};
for (i=0; i<nneg; i++) {
if (i%10==0) cerr<<nneg-i<<'\r';
AliITStrackV2 *ntrk=(AliITStrackV2 *)negtrks.UncheckedAt(i);
+
for (Int_t k=0; k<npos; k++) {
AliITStrackV2 *ptrk=(AliITStrackV2 *)postrks.UncheckedAt(k);
// This function returns the DCA between two tracks
// The tracks will be moved to the point of DCA !
//--------------------------------------------------------------------
+ Double_t dy2=n->GetSigmaY2() + 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);
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;
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);
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 (dd<dm) break;
dt1*=0.5; dt2*=0.5;
if (div>512) {
return 1.e+33;
}
- return TMath::Sqrt(dm);
+ //return TMath::Sqrt(dm);
+ return TMath::Sqrt(dx*dx + dy*dy + dz*dz);
}
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
timer.Stop(); timer.Print();
} // event loop
- //delete gAlice;
- //gAlice=0;
+ delete gAlice; gAlice=0;
file->Close();
}
#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
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
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff