Updated code for tracking V2 (from Y. Belikov)
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 29 Jan 2002 07:10:44 +0000 (07:10 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 29 Jan 2002 07:10:44 +0000 (07:10 +0000)
20 files changed:
ITS/AliCascadeComparison.C [new file with mode: 0644]
ITS/AliCascadeFindVertices.C [new file with mode: 0644]
ITS/AliCascadeVertex.cxx [new file with mode: 0644]
ITS/AliCascadeVertex.h [new file with mode: 0644]
ITS/AliCascadeVertexer.cxx [new file with mode: 0644]
ITS/AliCascadeVertexer.h [new file with mode: 0644]
ITS/AliITSComparisonV2.C
ITS/AliITSFindClustersV2.C
ITS/AliITSrecoV2.h
ITS/AliITStestV2.C
ITS/AliITStrackV2.cxx
ITS/AliITStrackV2.h
ITS/AliITStrackerV2.cxx
ITS/AliITStrackerV2.h
ITS/AliV0Comparison.C
ITS/AliV0FindVertices.C
ITS/AliV0vertexer.cxx
ITS/ITSHitsToFastPoints.C
ITS/ITSLinkDef.h
ITS/Makefile

diff --git a/ITS/AliCascadeComparison.C b/ITS/AliCascadeComparison.C
new file mode 100644 (file)
index 0000000..92758ff
--- /dev/null
@@ -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 <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;
+}
diff --git a/ITS/AliCascadeFindVertices.C b/ITS/AliCascadeFindVertices.C
new file mode 100644 (file)
index 0000000..44eb584
--- /dev/null
@@ -0,0 +1,42 @@
+#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;
+}
diff --git a/ITS/AliCascadeVertex.cxx b/ITS/AliCascadeVertex.cxx
new file mode 100644 (file)
index 0000000..d732940
--- /dev/null
@@ -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 <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;
+}
diff --git a/ITS/AliCascadeVertex.h b/ITS/AliCascadeVertex.h
new file mode 100644 (file)
index 0000000..6979e73
--- /dev/null
@@ -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 <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
+
+
diff --git a/ITS/AliCascadeVertexer.cxx b/ITS/AliCascadeVertexer.cxx
new file mode 100644 (file)
index 0000000..c692598
--- /dev/null
@@ -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 <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;
+}
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ITS/AliCascadeVertexer.h b/ITS/AliCascadeVertexer.h
new file mode 100644 (file)
index 0000000..3dae6f3
--- /dev/null
@@ -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
+
+
index 2f82cc2..04fe342 100644 (file)
@@ -133,7 +133,7 @@ Int_t AliITSComparisonV2(Int_t event=0) {
        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);
index d915e30..2842449 100644 (file)
   #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);
 
@@ -117,7 +57,7 @@ Int_t AliITSFindClustersV2(char SlowOrFast='s') {
 
    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();
@@ -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 ="<<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);
        }
index ef11775..16b91ce 100644 (file)
@@ -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;
index 91160ce..afe9257 100644 (file)
@@ -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;
 
index 4c36090..80f9a16 100644 (file)
@@ -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 <TMatrixD.h>
@@ -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; 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 {
@@ -142,12 +126,12 @@ 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;
 }
 
@@ -159,7 +143,7 @@ GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
   //------------------------------------------------------------------
   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";
@@ -276,14 +260,44 @@ Double_t x0) const {
 }
 
 //____________________________________________________________________________
-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";
@@ -338,29 +352,7 @@ AliITStrackV2::PropagateTo(Double_t xk,Double_t x0,Double_t rho) {
 
   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;
 }
@@ -435,124 +427,41 @@ Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
   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
   //------------------------------------------------------------------
@@ -570,7 +479,7 @@ AliITStrackV2::Propagate(Double_t alp,Double_t xk,Double_t x0,Double_t rho) {
   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";
@@ -589,11 +498,11 @@ AliITStrackV2::Propagate(Double_t alp,Double_t xk,Double_t x0,Double_t rho) {
 
   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";
@@ -672,28 +581,6 @@ AliITStrackV2::Propagate(Double_t alp,Double_t xk,Double_t x0,Double_t rho) {
   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;
@@ -717,9 +604,48 @@ Double_t AliITStrackV2::GetD() const {
   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  
@@ -729,7 +655,8 @@ Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
   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());
@@ -758,75 +685,7 @@ Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
   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 :)  
@@ -842,25 +701,15 @@ void AliITStrackV2::ResetCovariance() {
 
 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++) {
@@ -871,22 +720,12 @@ void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
     }
   } 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);
 }
-  
-
index 012c961..608b438 100644 (file)
@@ -7,6 +7,7 @@
 //                       ITS Track Class
 //
 //        Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
+//     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
 //-------------------------------------------------------------------------
 
 
@@ -36,16 +37,15 @@ public:
   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); }
@@ -95,7 +95,7 @@ private:
 
   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
 };
@@ -108,21 +108,19 @@ void AliITStrackV2::GetExternalParameters(Double_t& xr, Double_t x[5]) const {
      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
 
 
index d161efd..55e84ee 100644 (file)
@@ -17,6 +17,7 @@
 //               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>
@@ -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; 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;
@@ -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<<tpcLabel<<" *****************\n";
 #endif
+
        try {
          ResetTrackToFollow(*t);
        } catch (const Char_t *msg) {
@@ -211,43 +225,6 @@ cout<<tpcLabel<<" *****************\n";
        }
        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();
        }
@@ -355,24 +332,25 @@ for (Int_t k=0; k<nc; k++) {
          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) {
@@ -395,7 +373,8 @@ for (Int_t k=0; k<nc; k++) {
 
            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;
@@ -406,41 +385,18 @@ for (Int_t k=0; k<nc; k++) {
          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()); 
@@ -465,6 +421,7 @@ for (Int_t k=0; k<nc; k++) {
   return 0;
 }
 
+
 AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
   //--------------------------------------------------------------------
   //       Return pointer to a given cluster
@@ -482,19 +439,22 @@ void AliITStrackerV2::FollowProlongation() {
   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
@@ -502,10 +462,9 @@ cout<<fI<<' ';
     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: "
@@ -516,32 +475,36 @@ cout<<fI<<' ';
     //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;
@@ -562,20 +525,16 @@ cout<<fI<<' ';
      }
   }
 
-  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]);
@@ -583,22 +542,18 @@ Int_t AliITStrackerV2::TakeNextProlongation() {
   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";
@@ -609,11 +564,7 @@ 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
@@ -625,29 +576,34 @@ cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
 
   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
@@ -682,8 +638,9 @@ Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
   //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;}
@@ -769,40 +726,83 @@ cout<<np<<' '<<nz<<endl;
 }
 
 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);
 }
@@ -817,7 +817,7 @@ Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
   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();
 
index dc7f52e..b6ff039 100644 (file)
@@ -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() {
index 53bc2f3..92bdc39 100644 (file)
@@ -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 <iostream.h>
   #include <fstream.h>
@@ -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 (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();
index fa10eb1..6d1cf5c 100644 (file)
@@ -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
                    };
index 077fb22..435d25e 100644 (file)
@@ -78,6 +78,7 @@ Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) {
    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);
 
@@ -171,6 +172,12 @@ Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) {
   // 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);
@@ -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 (dd<dm) break;
         dt1*=0.5; dt2*=0.5;
         if (div>512) {
@@ -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);
 }
 
 
index c91441f..11d07d0 100644 (file)
@@ -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();
 }
 
index 2ee407b..4771f98 100644 (file)
 #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
index 2eeefa3..fd07018 100644 (file)
@@ -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