* provided "as is" without express or implied warranty. *
**************************************************************************/
+/* $Id$ */
+
//-------------------------------------------------------------------------
// Implementation of the ITS tracker class
// It reads AliITSclusterV2 clusters and creates AliITStrackMI tracks
// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
//
//-------------------------------------------------------------------------
-#include "AliITSrecoV2.h"
+
+#include <TMatrixD.h>
#include <TTree.h>
-#include "AliITSgeom.h"
+#include <TTreeStream.h>
+#include <TTree.h>
+
#include "AliESD.h"
+#include "AliESDV0MI.h"
+#include "AliHelix.h"
#include "AliITSclusterV2.h"
+#include "AliITSgeom.h"
#include "AliITStrackerMI.h"
-#include "TMatrixD.h"
-#include "TFile.h"
-#include "TTree.h"
-#include "AliHelix.h"
-#include "AliESDV0MI.h"
-
+#include "AliTrackPointArray.h"
+#include "AliAlignObj.h"
ClassImp(AliITStrackerMI)
//--------------------------------------------------------------------
//This is the AliITStrackerMI constructor
//--------------------------------------------------------------------
+ fCoeficients = 0;
+ fAfterV0 = kFALSE;
AliITSgeom *g=(AliITSgeom*)geom;
-
Float_t x,y,z;
Int_t i;
for (i=1; i<kMaxLayer+1; i++) {
for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
fLastLayerToTrackTo=kLastLayerToTrackTo;
+ for (Int_t i=0;i<100000;i++){
+ fBestTrackIndex[i]=0;
+ }
+ //
+ fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
+
+}
+AliITStrackerMI::~AliITStrackerMI()
+{
+ //
+ //destructor
+ //
+ if (fCoeficients) delete []fCoeficients;
+ if (fDebugStreamer) {
+ //fDebugStreamer->Close();
+ delete fDebugStreamer;
+ }
}
void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
// The clusters must be already loaded !
//--------------------------------------------------------------------
TObjArray itsTracks(15000);
-
+ fOriginal.Clear();
+ fEsd = event; // store pointer to the esd
{/* Read ESD tracks */
Int_t nentr=event->GetNumberOfTracks();
Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
try {
t=new AliITStrackMI(*esd);
} catch (const Char_t *msg) {
- Warning("Clusters2Tracks",msg);
+ //Warning("Clusters2Tracks",msg);
delete t;
continue;
}
// write expected q
t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.);
- if (TMath::Abs(t->fD[0])>10) {
- delete t;
- continue;
- }
-
- if (TMath::Abs(vdist)>20) {
- delete t;
- continue;
- }
- if (TMath::Abs(1/t->Get1Pt())<0.120) {
- delete t;
- continue;
+ if (esd->GetV0Index(0)>0 && t->fD[0]<30){
+ //track - can be V0 according to TPC
}
-
- if (CorrectForDeadZoneMaterial(t)!=0) {
- Warning("Clusters2Tracks",
- "failed to correct for the material in the dead zone !\n");
- delete t;
- continue;
+ else{
+ if (TMath::Abs(t->fD[0])>10) {
+ delete t;
+ continue;
+ }
+
+ if (TMath::Abs(vdist)>20) {
+ delete t;
+ continue;
+ }
+ if (TMath::Abs(1/t->Get1Pt())<0.120) {
+ delete t;
+ continue;
+ }
+
+ if (CorrectForDeadZoneMaterial(t)!=0) {
+ //Warning("Clusters2Tracks",
+ // "failed to correct for the material in the dead zone !\n");
+ delete t;
+ continue;
+ }
}
t->fReconstructed = kFALSE;
itsTracks.AddLast(t);
+ fOriginal.AddLast(t);
}
} /* End Read ESD tracks */
itsTracks.Sort();
+ fOriginal.Sort();
Int_t nentr=itsTracks.GetEntriesFast();
fTrackHypothesys.Expand(nentr);
+ fBestHypothesys.Expand(nentr);
MakeCoeficients(nentr);
Int_t ntrk=0;
for (fPass=0; fPass<2; fPass++) {
fI = 6;
ResetTrackToFollow(*t);
ResetBestTrack();
-
- FollowProlongationTree(t,i);
-
+ FollowProlongationTree(t,i,fConstraint[fPass]);
SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
//
CookdEdx(besttrack);
besttrack->fFakeRatio=1.;
CookLabel(besttrack,0.); //For comparison only
- // besttrack->UpdateESDtrack(AliESDtrack::kITSin);
- //
-
UpdateESDtrack(besttrack,AliESDtrack::kITSin);
-
- if ( besttrack->GetNumberOfClusters()<5 && fConstraint[fPass]) {
+
+ /*
+ if ( besttrack->GetNumberOfClusters()<6 && fConstraint[fPass]) {
continue;
}
if (besttrack->fChi2MIP[0]+besttrack->fNUsed>3.5) continue;
if ( (TMath::Abs(besttrack->fD[0]*besttrack->fD[0]+besttrack->fD[1]*besttrack->fD[1])>0.1) && fConstraint[fPass]) continue;
//delete itsTracks.RemoveAt(i);
+ */
+ if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
+
+
t->fReconstructed = kTRUE;
ntrk++;
}
}
//GetBestHypothesysMIP(itsTracks);
- FindV0(event);
-
+ UpdateTPCV0(event);
+ FindV02(event);
+ fAfterV0 = kTRUE;
+ //GetBestHypothesysMIP(itsTracks);
+ //
itsTracks.Delete();
//
Int_t entries = fTrackHypothesys.GetEntriesFast();
}
fTrackHypothesys.Delete();
+ fBestHypothesys.Delete();
+ fOriginal.Clear();
+ delete []fCoeficients;
+ fCoeficients=0;
Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
-
+
return 0;
}
try {
t=new AliITStrackMI(*esd);
} catch (const Char_t *msg) {
- Warning("PropagateBack",msg);
+ //Warning("PropagateBack",msg);
delete t;
continue;
}
fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
if (RefitAt(49.,&fTrackToFollow,t)) {
if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
- Warning("PropagateBack",
- "failed to correct for the material in the dead zone !\n");
+ //Warning("PropagateBack",
+ // "failed to correct for the material in the dead zone !\n");
delete t;
continue;
}
// "inward propagated" TPC tracks
// The clusters must be loaded !
//--------------------------------------------------------------------
+ RefitV02(event);
Int_t nentr=event->GetNumberOfTracks();
Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
try {
t=new AliITStrackMI(*esd);
} catch (const Char_t *msg) {
- Warning("RefitInward",msg);
+ //Warning("RefitInward",msg);
delete t;
continue;
}
t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.);
if (CorrectForDeadZoneMaterial(t)!=0) {
- Warning("RefitInward",
- "failed to correct for the material in the dead zone !\n");
+ //Warning("RefitInward",
+ // "failed to correct for the material in the dead zone !\n");
delete t;
continue;
}
CookLabel(&fTrackToFollow,0.0); //For comparison only
if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
- Double_t a=fTrackToFollow.GetAlpha();
- Double_t cs=TMath::Cos(a),sn=TMath::Sin(a);
- Double_t xv= GetX()*cs + GetY()*sn;
- Double_t yv=-GetX()*sn + GetY()*cs;
-
- Double_t c=fTrackToFollow.GetC(), snp=fTrackToFollow.GetSnp();
- Double_t x=fTrackToFollow.GetX(), y=fTrackToFollow.GetY();
- Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp));
- Double_t fv=TMath::ATan(tgfv);
-
- cs=TMath::Cos(fv); sn=TMath::Sin(fv);
- x = xv*cs + yv*sn;
- yv=-xv*sn + yv*cs; xv=x;
-
- if (fTrackToFollow.Propagate(fv+a,xv)) {
- fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
- //UseClusters(&fTrackToFollow);
- {
- AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ());
- c.SetSigmaY2(GetSigmaY()*GetSigmaY());
- c.SetSigmaZ2(GetSigmaZ()*GetSigmaZ());
- Double_t chi2=fTrackToFollow.GetPredictedChi2(&c);
- //Double_t chi2=GetPredictedChi2MI(&fTrackToFollow,&c,fI);
- if (chi2<kMaxChi2)
- if (fTrackToFollow.Update(&c,-chi2,0))
- //if (UpdateMI(&fTrackToFollow,&c,-chi2,0))
- fTrackToFollow.SetConstrainedESDtrack(chi2);
- }
- ntrk++;
- }
+ AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
+ esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit);
+ Float_t r[3]={0.,0.,0.};
+ Double_t maxD=3.;
+ esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
+ ntrk++;
}
}
delete t;
return fgLayers[l].GetCluster(c);
}
+Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
+ //
+ // Get track space point with index i
+ //
+ Int_t l=(index & 0xf0000000) >> 28;
+ Int_t c=(index & 0x0fffffff) >> 00;
+ AliITSclusterV2 *cl = fgLayers[l].GetCluster(c);
+ Int_t idet = cl->GetDetectorIndex();
+ const AliITSdetector &det = fgLayers[l].GetDetector(idet);
+ Float_t phi = det.GetPhi();
+ Float_t r = det.GetR();
+ Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
+ Float_t xyz[3];
+ xyz[0] = r*cp - cl->GetY()*sp;
+ xyz[1] = r*sp + cl->GetY()*cp;
+ xyz[2] = cl->GetZ();
+ p.SetXYZ(xyz[0],xyz[1],xyz[2]);
+ AliAlignObj::ELayerID iLayer = AliAlignObj::kInvalidLayer;
+ switch (l) {
+ case 0:
+ iLayer = AliAlignObj::kSPD1;
+ break;
+ case 1:
+ iLayer = AliAlignObj::kSPD2;
+ break;
+ case 2:
+ iLayer = AliAlignObj::kSDD1;
+ break;
+ case 3:
+ iLayer = AliAlignObj::kSDD2;
+ break;
+ case 4:
+ iLayer = AliAlignObj::kSSD1;
+ break;
+ case 5:
+ iLayer = AliAlignObj::kSSD2;
+ break;
+ default:
+ AliWarning(Form("Wrong layer index in ITS (%d) !",l));
+ break;
+ };
+ UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,idet);
+ p.SetVolumeID((UShort_t)volid);
+ return kTRUE;
+}
-void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex)
+void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
{
//--------------------------------------------------------------------
// Follow prolongation tree
//--------------------------------------------------------------------
+ //
+ AliESDtrack * esd = otrack->fESDtrack;
+ if (esd->GetV0Index(0)>0){
+ //
+ // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
+ // mapping of esd track is different as its track in Containers
+ // Need something more stable
+ // Indexes are set back againg to the ESD track indexes in UpdateTPCV0
+ for (Int_t i=0;i<3;i++){
+ Int_t index = esd->GetV0Index(i);
+ if (index==0) break;
+ AliESDV0MI * vertex = fEsd->GetV0MI(index);
+ if (vertex->GetStatus()<0) continue; // rejected V0
+ //
+ if (esd->GetSign()>0) {
+ vertex->SetIndex(0,esdindex);
+ }
+ else{
+ vertex->SetIndex(1,esdindex);
+ }
+ }
+ }
+ TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
+ if (!bestarray){
+ bestarray = new TObjArray(5);
+ fBestHypothesys.AddAt(bestarray,esdindex);
+ }
+ //
//setup tree of the prolongations
//
static AliITStrackMI tracks[7][100];
otrack->fNSkipped=0;
new (&(tracks[6][0])) AliITStrackMI(*otrack);
ntracks[6]=1;
- nindexes[6][0]=0;
+ for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
//
//
// follow prolongations
//
Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]));
Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]));
- if (fConstraint[fPass]){
+ if (constrain){
msy/=60; msz/=60.;
}
else{
Double_t x0;
Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0);
updatetrack->CorrectForMaterial(d,x0);
- if (fConstraint[fPass]) {
- updatetrack->fConstrain = fConstraint[fPass];
+ if (constrain) {
+ updatetrack->fConstrain = constrain;
fI = ilayer;
Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
Double_t xyz[]={GetX(),GetY(),GetZ()};
ntracks[ilayer]++;
} // create new hypothesy
} // loop over possible cluster prolongation
- // if (fConstraint[fPass]&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0){
- if (itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){
+ // if (constrain&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0){
+ if (constrain&&itrack<2&¤ttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){
AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
vtrack->fClIndex[ilayer]=0;
fI = ilayer;
vtrack->fNSkipped++;
ntracks[ilayer]++;
}
+
+ if (constrain&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){ //big theta -- for low mult. runs
+ AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
+ vtrack->fClIndex[ilayer]=0;
+ fI = ilayer;
+ Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
+ Double_t xyz[]={GetX(),GetY(),GetZ()};
+ Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
+ vtrack->Improve(d,xyz,ers);
+ vtrack->fNDeadZone++;
+ ntracks[ilayer]++;
+ }
+
} //loop over track candidates
//
if ( normalizedchi2[itrack]<3+0.5*ilayer) golds++;
if (ilayer>4) accepted++;
else{
- if ( fConstraint[fPass] && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
- if (!fConstraint[fPass] && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
+ if ( constrain && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
+ if (!constrain && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
}
}
TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
if (ntracks[ilayer]>90) ntracks[ilayer]=90;
} //loop over layers
//printf("%d\t%d\t%d\t%d\t%d\t%d\n",ntracks[0],ntracks[1],ntracks[2],ntracks[3],ntracks[4],ntracks[5]);
+ Int_t max = constrain? 20: 5;
- for (Int_t i=0;i<TMath::Min(20,ntracks[0]);i++) {
+ for (Int_t i=0;i<TMath::Min(max,ntracks[0]);i++) {
AliITStrackMI & track= tracks[0][nindexes[0][i]];
- if (!fConstraint[fPass]&&track.fNormChi2[0]>7.)continue;
+ if (track.GetNumberOfClusters()<2) continue;
+ if (!constrain&&track.fNormChi2[0]>7.)continue;
AddTrackHypothesys(new AliITStrackMI(track), esdindex);
}
- for (Int_t i=0;i<TMath::Min(4,ntracks[1]);i++) {
+ for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
AliITStrackMI & track= tracks[1][nindexes[1][i]];
if (track.GetNumberOfClusters()<4) continue;
- if (!fConstraint[fPass]&&track.fNormChi2[1]>7.)continue;
- if (fConstraint[fPass]) track.fNSkipped+=1;
- if (!fConstraint[fPass]) {
+ if (!constrain&&track.fNormChi2[1]>7.)continue;
+ if (constrain) track.fNSkipped+=1;
+ if (!constrain) {
track.fD[0] = track.GetD(GetX(),GetY());
track.fNSkipped+=4./(4.+8.*TMath::Abs(track.fD[0]));
if (track.fN+track.fNDeadZone+track.fNSkipped>6) {
}
//}
- if (!fConstraint[fPass]){
- for (Int_t i=0;i<TMath::Min(3,ntracks[2]);i++) {
+ if (!constrain){
+ for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
AliITStrackMI & track= tracks[2][nindexes[2][i]];
- if (track.GetNumberOfClusters()<4) continue;
- if (!fConstraint[fPass]&&track.fNormChi2[2]>7.)continue;
- if (fConstraint[fPass]) track.fNSkipped+=2;
- if (!fConstraint[fPass]){
+ if (track.GetNumberOfClusters()<3) continue;
+ if (!constrain&&track.fNormChi2[2]>7.)continue;
+ if (constrain) track.fNSkipped+=2;
+ if (!constrain){
track.fD[0] = track.GetD(GetX(),GetY());
track.fNSkipped+= 7./(7.+8.*TMath::Abs(track.fD[0]));
if (track.fN+track.fNDeadZone+track.fNSkipped>6) {
AddTrackHypothesys(new AliITStrackMI(track), esdindex);
}
}
+
+ if (!constrain){
+ //
+ // register best tracks - important for V0 finder
+ //
+ for (Int_t ilayer=0;ilayer<5;ilayer++){
+ if (ntracks[ilayer]==0) continue;
+ AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
+ if (track.GetNumberOfClusters()<1) continue;
+ CookLabel(&track,0);
+ bestarray->AddAt(new AliITStrackMI(track),ilayer);
+ }
+ }
+ //
+ // update TPC V0 information
+ //
+ if (otrack->fESDtrack->GetV0Index(0)>0){
+ Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
+ for (Int_t i=0;i<3;i++){
+ Int_t index = otrack->fESDtrack->GetV0Index(i);
+ if (index==0) break;
+ AliESDV0MI * vertex = fEsd->GetV0MI(index);
+ if (vertex->GetStatus()<0) continue; // rejected V0
+ //
+ if (otrack->fP4>0) {
+ vertex->SetIndex(0,esdindex);
+ }
+ else{
+ vertex->SetIndex(1,esdindex);
+ }
+ //find nearest layer with track info
+ Int_t nearestold = GetNearestLayer(vertex->GetXrp());
+ Int_t nearest = nearestold;
+ for (Int_t ilayer =nearest;ilayer<8;ilayer++){
+ if (ntracks[nearest]==0){
+ nearest = ilayer;
+ }
+ }
+ //
+ AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
+ if (nearestold<5&&nearest<5){
+ Bool_t accept = track.fNormChi2[nearest]<10;
+ if (accept){
+ if (track.fP4>0) {
+ vertex->SetP(track);
+ vertex->Update(fprimvertex);
+ // vertex->SetIndex(0,track.fESDtrack->GetID());
+ if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
+ }else{
+ vertex->SetM(track);
+ vertex->Update(fprimvertex);
+ //vertex->SetIndex(1,track.fESDtrack->GetID());
+ if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
+ }
+ vertex->SetStatus(vertex->GetStatus()+1);
+ }else{
+ // vertex->SetStatus(-2); // reject V0 - not enough clusters
+ }
+ }
+ // if (nearestold>3){
+// Int_t indexlayer = (ntracks[0]>0)? 0:1;
+// if (ntracks[indexlayer]>0){
+// AliITStrackMI & track= tracks[indexlayer][nindexes[indexlayer][0]];
+// if (track.GetNumberOfClusters()>4&&track.fNormChi2[indexlayer]<4){
+// vertex->SetStatus(-1); // reject V0 - clusters before
+// }
+// }
+// }
+ }
+ }
}
if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
}
+
Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
//--------------------------------------------------------------------
//This function adds a cluster to this layer
return 1;
}
fCurrentSlice=-1;
- if (fN==0) {fClusters[fN++]=c; return 0;}
- Int_t i=FindClusterIndex(c->GetZ());
- memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
- memmove(fY+i+1 ,fY+i,(fN-i)*sizeof(Float_t));
- memmove(fZ+i+1 ,fZ+i,(fN-i)*sizeof(Float_t));
- fClusters[i]=c; fN++;
- //
+ fClusters[fN]=c;
+ fN++;
AliITSdetector &det=GetDetector(c->GetDetectorIndex());
- Double_t y=fR*det.GetPhi() + c->GetY();
- if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
- fY[i] = y;
- fZ[i] = c->GetZ();
if (c->GetY()<det.GetYmin()) det.SetYmin(c->GetY());
if (c->GetY()>det.GetYmax()) det.SetYmax(c->GetY());
if (c->GetZ()<det.GetZmin()) det.SetZmin(c->GetZ());
//
//sort clusters
//
+ AliITSclusterV2 **clusters = new AliITSclusterV2*[fN];
+ Float_t *z = new Float_t[fN];
+ Int_t * index = new Int_t[fN];
+ //
+ for (Int_t i=0;i<fN;i++){
+ z[i] = fClusters[i]->GetZ();
+ }
+ TMath::Sort(fN,z,index,kFALSE);
+ for (Int_t i=0;i<fN;i++){
+ clusters[i] = fClusters[index[i]];
+ }
+ //
+ for (Int_t i=0;i<fN;i++){
+ fClusters[i] = clusters[i];
+ fZ[i] = fClusters[i]->GetZ();
+ AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
+ Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
+ if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
+ fY[i] = y;
+ }
+ delete[] index;
+ delete[] z;
+ delete[] clusters;
+ //
+
fYB[0]=10000000;
fYB[1]=-10000000;
for (Int_t i=0;i<fN;i++){
for (Int_t i=0;i<21;i++) {fBy20[i][0] = fYB[0]+(i-0.75)*fDy20; fBy20[i][1] = fYB[0]+(i+0.75)*fDy20;}
//
//
- for (Int_t i=0;i<fN;i++) for (Int_t irot=-1;irot<=1;irot++){
- Float_t curY = fY[i]+irot*2.*TMath::Pi()*fR;
- if (curY<fYB[0]-fDy5) continue;
- if (curY>fYB[1]+fDy5) continue;
- //
- // slice 10
- if (TMath::Abs(fYB[1]-fYB[0])<=0) continue;
- Float_t fslice = TMath::Nint(10.*(curY-fYB[0])/(fYB[1]-fYB[0]));
- Float_t ymiddle = fYB[0]+fslice*fDy10;
- for (Int_t di =-1;di<=1;di++){
- if (TMath::Abs(curY-(ymiddle+(float)di*fDy10))<0.75*fDy10){
- //
- Int_t slice = int(fslice+21.0001)-21+di;
- if (slice<0) continue;
- if (slice>10) continue;
- if (fN10[slice]>=kMaxClusterPerLayer10) break;
- fClusters10[slice][fN10[slice]] = fClusters[i];
- fY10[slice][fN10[slice]] = curY;
- fZ10[slice][fN10[slice]] = fZ[i];
- fClusterIndex10[slice][fN10[slice]]=i;
- fN10[slice]++;
+ for (Int_t i=0;i<fN;i++)
+ for (Int_t irot=-1;irot<=1;irot++){
+ Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
+ // slice 5
+ for (Int_t slice=0; slice<6;slice++){
+ if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<kMaxClusterPerLayer5){
+ fClusters5[slice][fN5[slice]] = fClusters[i];
+ fY5[slice][fN5[slice]] = curY;
+ fZ5[slice][fN5[slice]] = fZ[i];
+ fClusterIndex5[slice][fN5[slice]]=i;
+ fN5[slice]++;
+ }
}
- }
- //
- // slice 5
- fslice = TMath::Nint(5.*(curY-fYB[0])/(fYB[1]-fYB[0]));
- ymiddle = fYB[0]+fslice*fDy5;
- for (Int_t di =-1;di<=1;di++){
- if (TMath::Abs(curY-(ymiddle+(float)di*fDy5))<0.75*fDy5){
- //
- Int_t slice = int(fslice+21.0001)-21+di;
- if (slice<0) continue;
- if (slice>5) continue;
- if (fN5[slice]>=kMaxClusterPerLayer5) break;
- fClusters5[slice][fN5[slice]] = fClusters[i];
- fY5[slice][fN5[slice]] = curY;
- fZ5[slice][fN5[slice]] = fZ[i];
- fClusterIndex5[slice][fN5[slice]]=i;
- fN5[slice]++;
+ // slice 10
+ for (Int_t slice=0; slice<11;slice++){
+ if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<kMaxClusterPerLayer10){
+ fClusters10[slice][fN10[slice]] = fClusters[i];
+ fY10[slice][fN10[slice]] = curY;
+ fZ10[slice][fN10[slice]] = fZ[i];
+ fClusterIndex10[slice][fN10[slice]]=i;
+ fN10[slice]++;
+ }
}
+ // slice 20
+ for (Int_t slice=0; slice<21;slice++){
+ if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<kMaxClusterPerLayer20){
+ fClusters20[slice][fN20[slice]] = fClusters[i];
+ fY20[slice][fN20[slice]] = curY;
+ fZ20[slice][fN20[slice]] = fZ[i];
+ fClusterIndex20[slice][fN20[slice]]=i;
+ fN20[slice]++;
+ }
+ }
}
- //
- // slice 20
- fslice = TMath::Nint(20.*(curY-fYB[0])/(fYB[1]-fYB[0]));
- ymiddle = fYB[0]+fslice*fDy20;
- for (Int_t di =-1;di<=1;di++){
- if (TMath::Abs(curY-(ymiddle+(float)di*fDy20))<0.75*fDy20){
- //
- Int_t slice = int(fslice+21.0001)-21+di;
- if (slice<0) continue;
- if (slice>20) continue;
- if (fN20[slice]>=kMaxClusterPerLayer20) break;
- fClusters20[slice][fN20[slice]] = fClusters[i];
- fY20[slice][fN20[slice]] = curY;
- fZ20[slice][fN20[slice]] = fZ[i];
- fClusterIndex20[slice][fN20[slice]]=i;
- fN20[slice]++;
- }
+
+ //
+ // consistency check
+ //
+ for (Int_t i=0;i<fN-1;i++){
+ if (fZ[i]>fZ[i+1]){
+ printf("Bugg\n");
+ }
+ }
+ //
+ for (Int_t slice=0;slice<21;slice++)
+ for (Int_t i=0;i<fN20[slice]-1;i++){
+ if (fZ20[slice][i]>fZ20[slice][i+1]){
+ printf("Bugg\n");
}
}
-}
+}
Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
}
return m;
}
-/*
-void AliITStrackerMI::AliITSlayer::
-SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
- //--------------------------------------------------------------------
- // This function sets the "window"
- //--------------------------------------------------------------------
- fI=FindClusterIndex(zmin); fZmax=zmax;
- fImax = TMath::Min(FindClusterIndex(zmax)+1,fN);
- Double_t circle=2*TMath::Pi()*fR;
- if (ymax>circle) { ymax-=circle; ymin-=circle; }
- fYmin=ymin; fYmax=ymax;
- fSkip = 0;
-}
-*/
void AliITStrackerMI::AliITSlayer::
fAccepted =0;
}
-/*
-const AliITSclusterV2 *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){
+
+
+
+Int_t AliITStrackerMI::AliITSlayer::
+FindDetectorIndex(Double_t phi, Double_t z) const {
//--------------------------------------------------------------------
- // This function returns clusters within the "window"
+ //This function finds the detector crossed by the track
//--------------------------------------------------------------------
- const AliITSclusterV2 *cluster=0;
- for (Int_t i=fI; i<fN; i++) {
- const AliITSclusterV2 *c=fClusters[i];
- if (c->GetZ() > fZmax) break;
- // if (c->IsUsed()) continue;
- const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
- Double_t y=fR*det.GetPhi() + c->GetY();
+ Double_t dphi=-(phi-fPhiOffset);
+ if (dphi < 0) dphi += 2*TMath::Pi();
+ else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
+ Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
+ if (np>=fNladders) np-=fNladders;
+ if (np<0) np+=fNladders;
- if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
- if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
+ Double_t dz=fZOffset-z;
+ Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
+ if (nz>=fNdetectors) return -1;
+ if (nz<0) return -1;
- if (y<fYmin) continue;
- if (y>fYmax) continue;
- cluster=c; ci=i;
- fI=i+1;
- break;
- }
- return cluster;
+ return np*fNdetectors + nz;
}
-*/
+
const AliITSclusterV2 *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){
//--------------------------------------------------------------------
for (Int_t i=fI; i<fImax; i++) {
Double_t y = fY[i];
if (fYmax<y) y -= rpi2;
+ if (fYmin>y) y += rpi2;
if (y<fYmin) continue;
if (y>fYmax) continue;
if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
-Int_t AliITStrackerMI::AliITSlayer::
-FindDetectorIndex(Double_t phi, Double_t z) const {
- //--------------------------------------------------------------------
- //This function finds the detector crossed by the track
- //--------------------------------------------------------------------
- Double_t dphi=-(phi-fPhiOffset);
- if (dphi < 0) dphi += 2*TMath::Pi();
- else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
- Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
- if (np>=fNladders) np-=fNladders;
- if (np<0) np+=fNladders;
-
- Double_t dz=fZOffset-z;
- Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
- if (nz>=fNdetectors) return -1;
- if (nz<0) return -1;
-
- return np*fNdetectors + nz;
-}
-
Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
const {
//--------------------------------------------------------------------
Int_t idx=index[i];
if (idx>0) {
const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx);
- if (idet != c->GetDetectorIndex()) {
- idet=c->GetDetectorIndex();
- const AliITSdetector &det=layer.GetDetector(idet);
- if (!t->Propagate(det.GetPhi(),det.GetR())) {
- return kFALSE;
- }
- t->SetDetectorIndex(idet);
- }
- //Double_t chi2=t->GetPredictedChi2(c);
- Int_t layer = (idx & 0xf0000000) >> 28;;
- Double_t chi2=GetPredictedChi2MI(t,c,layer);
- if (chi2<maxchi2) {
- cl=c;
- maxchi2=chi2;
- } else {
- return kFALSE;
+ if (c){
+ if (idet != c->GetDetectorIndex()) {
+ idet=c->GetDetectorIndex();
+ const AliITSdetector &det=layer.GetDetector(idet);
+ if (!t->Propagate(det.GetPhi(),det.GetR())) {
+ return kFALSE;
+ }
+ t->SetDetectorIndex(idet);
+ }
+ //Double_t chi2=t->GetPredictedChi2(c);
+ Int_t layer = (idx & 0xf0000000) >> 28;;
+ Double_t chi2=GetPredictedChi2MI(t,c,layer);
+ if (chi2<maxchi2) {
+ cl=c;
+ maxchi2=chi2;
+ } else {
+ return kFALSE;
+ }
}
}
/*
}
-Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
-{
+Bool_t
+AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,const Int_t *clindex) {
+ //--------------------------------------------------------------------
+ // This function refits the track "t" at the position "x" using
+ // the clusters from array
+ //--------------------------------------------------------------------
+ Int_t index[kMaxLayer];
+ Int_t k;
+ for (k=0; k<kMaxLayer; k++) index[k]=-1;
//
- // calculate normalized chi2
- // return NormalizedChi2(track,0);
- Float_t chi2 = 0;
- Float_t sum=0;
- Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
- // track->fdEdxMismatch=0;
- Float_t dedxmismatch =0;
- Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
- if (mode<100){
- for (Int_t i = 0;i<6;i++){
- if (track->fClIndex[i]>0){
- Float_t cerry, cerrz;
- if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
- else
- { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
- cerry*=cerry;
- cerrz*=cerrz;
- Float_t cchi2 = (track->fDy[i]*track->fDy[i]/cerry)+(track->fDz[i]*track->fDz[i]/cerrz);
- if (i>1){
- Float_t ratio = track->fNormQ[i]/track->fExpQ;
- if (ratio<0.5) {
- cchi2+=(0.5-ratio)*10.;
- //track->fdEdxMismatch+=(0.5-ratio)*10.;
- dedxmismatch+=(0.5-ratio)*10.;
- }
- }
- if (i<2 ||i>3){
- AliITSclusterV2 * cl = (AliITSclusterV2*)GetCluster( track->fClIndex[i]);
- Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
- if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
- if (i<2) chi2+=2*cl->GetDeltaProbability();
- }
- chi2+=cchi2;
- sum++;
- }
- }
- if (TMath::Abs(track->fdEdxMismatch-dedxmismatch)>0.0001){
- track->fdEdxMismatch = dedxmismatch;
- }
+ for (k=0; k<kMaxLayer; k++) {
+ index[k]=clindex[k];
}
- else{
- for (Int_t i = 0;i<4;i++){
- if (track->fClIndex[i]>0){
- Float_t cerry, cerrz;
- if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
- else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
- cerry*=cerry;
- cerrz*=cerrz;
- chi2+= (track->fDy[i]*track->fDy[i]/cerry);
- chi2+= (track->fDz[i]*track->fDz[i]/cerrz);
- sum++;
- }
- }
- for (Int_t i = 4;i<6;i++){
- if (track->fClIndex[i]>0){
- Float_t cerry, cerrz;
+
+ Int_t from, to, step;
+ if (xx > t->GetX()) {
+ from=0; to=kMaxLayer;
+ step=+1;
+ } else {
+ from=kMaxLayer-1; to=-1;
+ step=-1;
+ }
+
+ for (Int_t i=from; i != to; i += step) {
+ AliITSlayer &layer=fgLayers[i];
+ Double_t r=layer.GetR();
+ if (step<0 && xx>r) break; //
+ {
+ Double_t hI=i-0.5*step;
+ if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
+ Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
+ Double_t d=0.0034, x0=38.6;
+ if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
+ if (!t->PropagateTo(rs,-step*d,x0)) {
+ return kFALSE;
+ }
+ }
+ }
+
+ // remember old position [SR, GSI 18.02.2003]
+ Double_t oldX=0., oldY=0., oldZ=0.;
+ if (t->IsStartedTimeIntegral() && step==1) {
+ t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
+ }
+ //
+
+ Double_t x,y,z;
+ if (!t->GetGlobalXYZat(r,x,y,z)) {
+ return kFALSE;
+ }
+ Double_t phi=TMath::ATan2(y,x);
+ Int_t idet=layer.FindDetectorIndex(phi,z);
+ if (idet<0) {
+ return kFALSE;
+ }
+ const AliITSdetector &det=layer.GetDetector(idet);
+ phi=det.GetPhi();
+ if (!t->Propagate(phi,det.GetR())) {
+ return kFALSE;
+ }
+ t->SetDetectorIndex(idet);
+
+ const AliITSclusterV2 *cl=0;
+ Double_t maxchi2=1000.*kMaxChi2;
+
+ Int_t idx=index[i];
+ if (idx>0) {
+ const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx);
+ if (c){
+ if (idet != c->GetDetectorIndex()) {
+ idet=c->GetDetectorIndex();
+ const AliITSdetector &det=layer.GetDetector(idet);
+ if (!t->Propagate(det.GetPhi(),det.GetR())) {
+ return kFALSE;
+ }
+ t->SetDetectorIndex(idet);
+ }
+ //Double_t chi2=t->GetPredictedChi2(c);
+ Int_t layer = (idx & 0xf0000000) >> 28;;
+ Double_t chi2=GetPredictedChi2MI(t,c,layer);
+ if (chi2<maxchi2) {
+ cl=c;
+ maxchi2=chi2;
+ } else {
+ return kFALSE;
+ }
+ }
+ }
+ /*
+ if (cl==0)
+ if (t->GetNumberOfClusters()>2) {
+ Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
+ Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
+ Double_t zmin=t->GetZ() - dz;
+ Double_t zmax=t->GetZ() + dz;
+ Double_t ymin=t->GetY() + phi*r - dy;
+ Double_t ymax=t->GetY() + phi*r + dy;
+ layer.SelectClusters(zmin,zmax,ymin,ymax);
+
+ const AliITSclusterV2 *c=0; Int_t ci=-1;
+ while ((c=layer.GetNextCluster(ci))!=0) {
+ if (idet != c->GetDetectorIndex()) continue;
+ Double_t chi2=t->GetPredictedChi2(c);
+ if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
+ }
+ }
+ */
+ if (cl) {
+ //if (!t->Update(cl,maxchi2,idx)) {
+ if (!UpdateMI(t,cl,maxchi2,idx)) {
+ return kFALSE;
+ }
+ t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
+ }
+
+ {
+ Double_t x0;
+ Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
+ t->CorrectForMaterial(-step*d,x0);
+ }
+
+ // track time update [SR, GSI 17.02.2003]
+ if (t->IsStartedTimeIntegral() && step==1) {
+ Double_t newX, newY, newZ;
+ t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
+ Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
+ (oldZ-newZ)*(oldZ-newZ);
+ t->AddTimeStep(TMath::Sqrt(dL2));
+ }
+ //
+
+ }
+
+ if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
+ return kTRUE;
+}
+
+
+
+Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
+{
+ //
+ // calculate normalized chi2
+ // return NormalizedChi2(track,0);
+ Float_t chi2 = 0;
+ Float_t sum=0;
+ Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
+ // track->fdEdxMismatch=0;
+ Float_t dedxmismatch =0;
+ Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
+ if (mode<100){
+ for (Int_t i = 0;i<6;i++){
+ if (track->fClIndex[i]>0){
+ Float_t cerry, cerrz;
+ if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
+ else
+ { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
+ cerry*=cerry;
+ cerrz*=cerrz;
+ Float_t cchi2 = (track->fDy[i]*track->fDy[i]/cerry)+(track->fDz[i]*track->fDz[i]/cerrz);
+ if (i>1){
+ Float_t ratio = track->fNormQ[i]/track->fExpQ;
+ if (ratio<0.5) {
+ cchi2+=(0.5-ratio)*10.;
+ //track->fdEdxMismatch+=(0.5-ratio)*10.;
+ dedxmismatch+=(0.5-ratio)*10.;
+ }
+ }
+ if (i<2 ||i>3){
+ AliITSclusterV2 * cl = (AliITSclusterV2*)GetCluster( track->fClIndex[i]);
+ Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
+ if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
+ if (i<2) chi2+=2*cl->GetDeltaProbability();
+ }
+ chi2+=cchi2;
+ sum++;
+ }
+ }
+ if (TMath::Abs(track->fdEdxMismatch-dedxmismatch)>0.0001){
+ track->fdEdxMismatch = dedxmismatch;
+ }
+ }
+ else{
+ for (Int_t i = 0;i<4;i++){
+ if (track->fClIndex[i]>0){
+ Float_t cerry, cerrz;
+ if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
+ else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
+ cerry*=cerry;
+ cerrz*=cerrz;
+ chi2+= (track->fDy[i]*track->fDy[i]/cerry);
+ chi2+= (track->fDz[i]*track->fDz[i]/cerrz);
+ sum++;
+ }
+ }
+ for (Int_t i = 4;i<6;i++){
+ if (track->fClIndex[i]>0){
+ Float_t cerry, cerrz;
if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
cerry*=cerry;
{
//---------------------------------------------
// register track to the list
+ //
+ if (track->fESDtrack->GetKinkIndex(0)!=0) return; //don't register kink tracks
+ //
+ //
for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
Int_t index = track->GetClusterIndex(icluster);
Int_t l=(index & 0xf0000000) >> 28;
}
}
else{
- delete array->RemoveAt(itrack);
+ if (track->fConstrain || track->fN>5){ //keep best short tracks - without vertex constrain
+ delete array->RemoveAt(itrack);
+ }
}
}
if (!besttrack) return;
track->fChi2MIP[0] = GetNormalizedChi2(track, mode);
if (track->fChi2MIP[0]<kMaxChi2PerCluster[0])
chi2[itrack] = track->fChi2MIP[0];
- else
- delete array->RemoveAt(itrack);
+ else{
+ if (track->fConstrain || track->fN>5){ //keep best short tracks - without vertex constrain
+ delete array->RemoveAt(itrack);
+ }
+ }
}
}
//
track->fChi2MIP[0] = GetNormalizedChi2(track,mode);
if (track->fChi2MIP[0]<kMaxChi2PerCluster[0])
chi2[itrack] = track->fChi2MIP[0]-0*(track->GetNumberOfClusters()+track->fNDeadZone);
- else
- delete array->RemoveAt(itrack);
+ else
+ {
+ if (track->fConstrain || track->fN>5){ //keep best short tracks - without vertex constrain
+ delete array->RemoveAt(itrack);
+ }
+ }
}
}
entries = array->GetEntriesFast();
//
+ //
if (entries>0){
TObjArray * newarray = new TObjArray();
TMath::Sort(entries,chi2,index,kFALSE);
if (!track) continue;
if (accepted>maxcut) break;
track->fChi2MIP[0] = GetNormalizedChi2(track,mode);
- if (track->GetNumberOfClusters()<6 && (track->fChi2MIP[0]+track->fNUsed>minchi2)){
- delete array->RemoveAt(index[i]);
- continue;
+ if (track->fConstrain || track->fN>5){ //keep best short tracks - without vertex constrain
+ if (track->GetNumberOfClusters()<6 && (track->fChi2MIP[0]+track->fNUsed>minchi2)){
+ delete array->RemoveAt(index[i]);
+ continue;
+ }
}
- if (track->fChi2MIP[0]+track->fNUsed<minchi2 && track->GetNumberOfClusters()>=minn){
- accepted++;
+ Bool_t shortbest = !track->fConstrain && track->fN<6;
+ if ((track->fChi2MIP[0]+track->fNUsed<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
+ if (!shortbest) accepted++;
//
newarray->AddLast(array->RemoveAt(index[i]));
for (Int_t i=0;i<6;i++){
//
AliITStrackMI * backtrack = new AliITStrackMI(*original);
AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
+ Double_t xyzv[]={GetX(),GetY(),GetZ()};
+ Double_t ersv[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
//
for (Int_t i=0;i<entries;i++){
AliITStrackMI * track = (AliITStrackMI*)array->At(i);
if (!track) continue;
+ Float_t sigmarfi,sigmaz;
+ GetDCASigma(track,sigmarfi,sigmaz);
+ track->fDnorm[0] = sigmarfi;
+ track->fDnorm[1] = sigmaz;
+ //
track->fChi2MIP[1] = 1000000;
track->fChi2MIP[2] = 1000000;
track->fChi2MIP[3] = 1000000;
//
// backtrack
backtrack = new(backtrack) AliITStrackMI(*track);
- backtrack->ResetCovariance();
- backtrack->ResetCovariance();
+ if (track->fConstrain){
+ if (!backtrack->PropagateTo(3.,0.0028,65.19)) continue;
+ if (!backtrack->Improve(0,xyzv,ersv)) continue;
+ if (!backtrack->PropagateTo(2.,0.0028,0)) continue;
+ if (!backtrack->Improve(0,xyzv,ersv)) continue;
+ if (!backtrack->PropagateTo(1.,0.0028,0)) continue;
+ if (!backtrack->Improve(0,xyzv,ersv)) continue;
+ if (!backtrack->PropagateToVertex()) continue;
+ backtrack->ResetCovariance();
+ if (!backtrack->Improve(0,xyzv,ersv)) continue;
+ }else{
+ backtrack->ResetCovariance();
+ }
backtrack->ResetClusters();
+
Double_t x = original->GetX();
if (!RefitAt(x,backtrack,track)) continue;
+ //
track->fChi2MIP[1] = NormalizedChi2(backtrack,0);
//for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
if (track->fChi2MIP[1]>kMaxChi2PerCluster[1]*6.) continue;
if (!(track->fConstrain)&&track->fChi2MIP[1]>kMaxChi2PerCluster[1]) continue;
Bool_t isOK=kTRUE;
- /*
- for (Int_t i=0;i<6;i++){
- if (track->fClIndex[i]>0){
- Double_t dy1 = (track->fDy[i]/track->fSigmaY[i]);
- Double_t dz1 = (track->fDz[i]/track->fSigmaZ[i]);
- Double_t dy2 = (backtrack->fDy[i]/backtrack->fSigmaY[i]);
- Double_t dz2 = (backtrack->fDz[i]/backtrack->fSigmaZ[i]);
- if (TMath::Min(dy1*dy1+dz1*dz1,dy2*dy2+dz2*dz2)> kMaxChi2sR[i]) isOK =kFALSE;
- track->fDy[i+6] = backtrack->fDy[i];
- track->fDz[i+6] = backtrack->fDz[i];
- track->fSigmaY[i+6] = backtrack->fSigmaY[i];
- track->fSigmaZ[i+6] = backtrack->fSigmaZ[i];
- }
- else{
- if (i==5){
- if (track->fClIndex[i-1]>0){
- Double_t dy2 = (backtrack->fDy[i-1]/backtrack->fSigmaY[i-1]);
- Double_t dz2 = (backtrack->fDz[i-1]/backtrack->fSigmaZ[i-1]);
- if (dy2*dy2+dz2*dz2> kMaxChi2sR[i]) isOK =kFALSE;
- }
- else isOK = kFALSE;
- }
- }
- }
- */
if(!isOK) continue;
//
//forward track - without constraint
forwardtrack = new(forwardtrack) AliITStrackMI(*original);
forwardtrack->ResetClusters();
x = track->GetX();
- if (!RefitAt(x,forwardtrack,track)) continue;
+ RefitAt(x,forwardtrack,track);
track->fChi2MIP[2] = NormalizedChi2(forwardtrack,0);
if (track->fChi2MIP[2]>kMaxChi2PerCluster[2]*6.0) continue;
if (!(track->fConstrain)&&track->fChi2MIP[2]>kMaxChi2PerCluster[2]) continue;
for (Int_t i=0;i<entries;i++){
AliITStrackMI * track = (AliITStrackMI*)array->At(i);
if (!track) continue;
+
if (accepted>checkmax || track->fChi2MIP[3]>kMaxChi2PerCluster[3]*6. ||
(track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
track->fChi2MIP[0]>besttrack->fChi2MIP[0]+2.*besttrack->fNUsed+3.){
- delete array->RemoveAt(i);
- continue;
+ if (track->fConstrain || track->fN>5){ //keep best short tracks - without vertex constrain
+ delete array->RemoveAt(i);
+ continue;
+ }
}
else{
accepted++;
array->Compress();
SortTrackHypothesys(esdindex,checkmax,1);
array = (TObjArray*) fTrackHypothesys.At(esdindex);
+ if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
besttrack = (AliITStrackMI*)array->At(0);
if (!besttrack) return 0;
besttrack->fChi2MIP[8]=0;
Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
if (besttrack->fConstrain&&besttrack->fChi2MIP[0]<kMaxChi2PerCluster[0]&&besttrack->fChi2MIP[1]<kMaxChi2PerCluster[1]
&&besttrack->fChi2MIP[2]<kMaxChi2PerCluster[2]&&besttrack->fChi2MIP[3]<kMaxChi2PerCluster[3]){
- //if (besttrack->fChi2MIP[0]<kMaxChi2PerCluster[0]&&besttrack->fChi2MIP[1]<kMaxChi2PerCluster[1]
- // &&besttrack->fChi2MIP[2]<kMaxChi2PerCluster[2]&&besttrack->fChi2MIP[3]<kMaxChi2PerCluster[3]){
RegisterClusterTracks(besttrack,esdindex);
}
//
if (shared>2.5) return 0;
if (shared>1.0) return besttrack;
//
- // don't sign clusters if not gold track
- Double_t deltad = besttrack->GetD(GetX(),GetY());
- Double_t deltaz = besttrack->GetZat(GetX()) - GetZ();
- Double_t deltaprim = TMath::Sqrt(deltad*deltad+deltaz*deltaz);
- if (deltaprim>0.1 && (fConstraint[fPass])) return besttrack;
- if (TMath::Abs(deltad)>0.1) return besttrack;
- if (TMath::Abs(deltaz)>0.1) return besttrack;
- if (besttrack->fChi2MIP[0]>4.) return besttrack;
- if (besttrack->fChi2MIP[1]>4.) return besttrack;
- if (besttrack->fChi2MIP[2]>4.) return besttrack;
- if (besttrack->fChi2MIP[3]>4.) return besttrack;
+ // Don't sign clusters if not gold track
+ //
+ if (!besttrack->IsGoldPrimary()) return besttrack;
+ if (besttrack->fESDtrack->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
+ //
if (fConstraint[fPass]){
//
// sign clusters
{
//
// get "best" hypothesys
- //for (Int_t ilayer=0;ilayer<6;ilayer++) fgLayers[ilayer].ResetWeights();
-
+ //
Int_t nentries = itsTracks.GetEntriesFast();
for (Int_t i=0;i<nentries;i++){
for (Int_t j=0;j<array->GetEntriesFast();j++){
AliITStrackMI* track = (AliITStrackMI*)array->At(j);
if (!track) continue;
+ if (track->fGoldV0) {
+ longtrack = track; //gold V0 track taken
+ break;
+ }
if (track->GetNumberOfClusters()+track->fNDeadZone<minn) continue;
- if (track->GetNumberOfClusters()+track->fNDeadZone>minn) maxchi2 = track->fChi2MIP[0];
- if (track->fChi2MIP[0]>maxchi2) continue;
+ Float_t chi2 = track->fChi2MIP[0];
+ if (fAfterV0){
+ if (!track->fGoldV0&&track->fConstrain==kFALSE) chi2+=5;
+ }
+ if (track->GetNumberOfClusters()+track->fNDeadZone>minn) maxchi2 = track->fChi2MIP[0];
+ //
+ if (chi2 > maxchi2) continue;
minn= track->GetNumberOfClusters()+track->fNDeadZone;
- maxchi2 = track->fChi2MIP[0];
+ maxchi2 = chi2;
longtrack=track;
- break;
}
+ //
+ //
+ //
AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
if (!longtrack) {longtrack = besttrack;}
else besttrack= longtrack;
+ //
if (besttrack){
Int_t list[6];
AliITSclusterV2 * clist[6];
besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
}
}
+ if (besttrack&&fAfterV0){
+ UpdateESDtrack(besttrack,AliESDtrack::kITSin);
+ }
if (besttrack&&fConstraint[fPass])
UpdateESDtrack(besttrack,AliESDtrack::kITSin);
//if (besttrack&&besttrack->fConstrain)
void AliITStrackerMI::MakeCoeficients(Int_t ntracks){
//
//
+ if (fCoeficients) delete []fCoeficients;
fCoeficients = new Float_t[ntracks*48];
for (Int_t i=0;i<ntracks*48;i++) fCoeficients[i]=-1.;
}
//STRIPS
if (layer>3){
+ //factor 1.8 appears in new simulation
+ //
+ Float_t scale=1.8;
if (cl->GetNy()==100||cl->GetNz()==100){
- erry = 0.004;
- errz = 0.2;
+ erry = 0.004*scale;
+ errz = 0.2*scale;
return 100;
}
if (cl->GetNy()+cl->GetNz()>12){
- erry = 0.06;
- errz = 0.57;
+ erry = 0.06*scale;
+ errz = 0.57*scale;
return 100;
}
Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
//
if (cl->GetType()==1 || cl->GetType()==10 ){
if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
- errz = 0.043;
- erry = 0.00094;
+ errz = 0.043*scale;
+ erry = 0.00094*scale;
return 101;
}
if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
- errz = 0.06;
- erry =0.0013;
+ errz = 0.06*scale;
+ erry =0.0013*scale;
return 102;
}
- erry = 0.0027;
- errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15);
+ erry = 0.0027*scale;
+ errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
return 103;
}
if (cl->GetType()==2 || cl->GetType()==11 ){
- erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05);
- errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5);
+ erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
+ errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
return 104;
}
if (cl->GetType()>100 ){
if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
- errz = 0.05;
- erry = 0.00096;
+ errz = 0.05*scale;
+ erry = 0.00096*scale;
return 105;
}
if (cl->GetNy()+cl->GetNz()-nz-ny<1){
- errz = 0.10;
- erry = 0.0025;
+ errz = 0.10*scale;
+ erry = 0.0025*scale;
return 106;
}
- errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4);
- erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05);
+ errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
+ erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
return 107;
}
Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
+void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
+{
+ //
+ //DCA sigmas parameterization
+ //to be paramterized using external parameters in future
+ //
+ //
+ sigmarfi = 0.004+1.4 *TMath::Abs(track->fP4)+332.*track->fP4*track->fP4;
+ sigmaz = 0.011+4.37*TMath::Abs(track->fP4);
+}
+
+
void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
{
//
AliITStrackMI * oldtrack = (AliITStrackMI*)(track->fESDtrack->GetITStrack());
if (oldtrack) delete oldtrack;
track->fESDtrack->SetITStrack(new AliITStrackMI(*track));
+ if (TMath::Abs(track->fDnorm[1])<0.000000001){
+ printf("Problem\n");
+ }
}
+Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
+ //
+ //Get nearest upper layer close to the point xr.
+ // rough approximation
+ //
+ const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
+ Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+ Int_t res =6;
+ for (Int_t i=0;i<6;i++){
+ if (radius<kRadiuses[i]){
+ res =i;
+ break;
+ }
+ }
+ return res;
+}
+void AliITStrackerMI::UpdateTPCV0(AliESD *event){
+ //
+ //try to update, or reject TPC V0s
+ //
+ Int_t nv0s = event->GetNumberOfV0MIs();
+ Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
+ for (Int_t i=0;i<nv0s;i++){
+ AliESDV0MI * vertex = event->GetV0MI(i);
+ Int_t ip = vertex->GetIndex(0);
+ Int_t im = vertex->GetIndex(1);
+ //
+ TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
+ TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
+ AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
+ AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
+ //
+ //
+ if (trackp){
+ if (trackp->fN+trackp->fNDeadZone>5.5){
+ if (trackp->fConstrain&&trackp->fChi2MIP[0]<3) vertex->SetStatus(-100);
+ if (!trackp->fConstrain&&trackp->fChi2MIP[0]<2) vertex->SetStatus(-100);
+ }
+ }
+
+ if (trackm){
+ if (trackm->fN+trackm->fNDeadZone>5.5){
+ if (trackm->fConstrain&&trackm->fChi2MIP[0]<3) vertex->SetStatus(-100);
+ if (!trackm->fConstrain&&trackm->fChi2MIP[0]<2) vertex->SetStatus(-100);
+ }
+ }
+ if (vertex->GetStatus()==-100) continue;
+ //
+ Int_t clayer = GetNearestLayer(vertex->GetXrp());
+ vertex->SetNBefore(clayer); //
+ vertex->SetChi2Before(9*clayer); //
+ vertex->SetNAfter(6-clayer); //
+ vertex->SetChi2After(0); //
+ //
+ if (clayer >1 ){ // calculate chi2 before vertex
+ Float_t chi2p = 0, chi2m=0;
+ //
+ if (trackp){
+ for (Int_t ilayer=0;ilayer<clayer;ilayer++){
+ if (trackp->fClIndex[ilayer]>0){
+ chi2p+=trackp->fDy[ilayer]*trackp->fDy[ilayer]/(trackp->fSigmaY[ilayer]*trackp->fSigmaY[ilayer])+
+ trackp->fDz[ilayer]*trackp->fDz[ilayer]/(trackp->fSigmaZ[ilayer]*trackp->fSigmaZ[ilayer]);
+ }
+ else{
+ chi2p+=9;
+ }
+ }
+ }else{
+ chi2p = 9*clayer;
+ }
+ //
+ if (trackm){
+ for (Int_t ilayer=0;ilayer<clayer;ilayer++){
+ if (trackm->fClIndex[ilayer]>0){
+ chi2m+=trackm->fDy[ilayer]*trackm->fDy[ilayer]/(trackm->fSigmaY[ilayer]*trackm->fSigmaY[ilayer])+
+ trackm->fDz[ilayer]*trackm->fDz[ilayer]/(trackm->fSigmaZ[ilayer]*trackm->fSigmaZ[ilayer]);
+ }
+ else{
+ chi2m+=9;
+ }
+ }
+ }else{
+ chi2m = 9*clayer;
+ }
+ vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
+ if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex
+ }
+
+ if (clayer < 5 ){ // calculate chi2 after vertex
+ Float_t chi2p = 0, chi2m=0;
+ //
+ if (trackp&&TMath::Abs(trackp->fP3)<1.){
+ for (Int_t ilayer=clayer;ilayer<6;ilayer++){
+ if (trackp->fClIndex[ilayer]>0){
+ chi2p+=trackp->fDy[ilayer]*trackp->fDy[ilayer]/(trackp->fSigmaY[ilayer]*trackp->fSigmaY[ilayer])+
+ trackp->fDz[ilayer]*trackp->fDz[ilayer]/(trackp->fSigmaZ[ilayer]*trackp->fSigmaZ[ilayer]);
+ }
+ else{
+ chi2p+=9;
+ }
+ }
+ }else{
+ chi2p = 0;
+ }
+ //
+ if (trackm&&TMath::Abs(trackm->fP3)<1.){
+ for (Int_t ilayer=clayer;ilayer<6;ilayer++){
+ if (trackm->fClIndex[ilayer]>0){
+ chi2m+=trackm->fDy[ilayer]*trackm->fDy[ilayer]/(trackm->fSigmaY[ilayer]*trackm->fSigmaY[ilayer])+
+ trackm->fDz[ilayer]*trackm->fDz[ilayer]/(trackm->fSigmaZ[ilayer]*trackm->fSigmaZ[ilayer]);
+ }
+ else{
+ chi2m+=9;
+ }
+ }
+ }else{
+ chi2m = 0;
+ }
+ vertex->SetChi2After(TMath::Max(chi2p,chi2m));
+ if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS
+ }
+ }
+ //
+}
-void AliITStrackerMI::FindV0(AliESD *event)
+
+void AliITStrackerMI::FindV02(AliESD *event)
{
//
- // fast V0 finder
+ // V0 finder
//
- //fV0Array->Clean();
- AliHelix helixes[30000];
- TObjArray trackarray(30000);
- Float_t dist[30000];
- Float_t mindist[30000];
- AliESDV0MI *vertexarray = new AliESDV0MI[20000];
- AliV0vertex *oldvertexarray = new AliV0vertex[20000];
- AliESDV0MI *pvertex = &vertexarray[0];
-
+ // Cuts on DCA - R dependent
+ // max distance DCA between 2 tracks cut
+ // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
+ //
+ const Float_t kMaxDist0 = 0.1;
+ const Float_t kMaxDist1 = 0.1;
+ const Float_t kMaxDist = 1;
+ const Float_t kMinPointAngle = 0.85;
+ const Float_t kMinPointAngle2 = 0.99;
+ const Float_t kMinR = 0.5;
+ const Float_t kMaxR = 220;
+ //const Float_t kCausality0Cut = 0.19;
+ //const Float_t kLikelihood01Cut = 0.25;
+ //const Float_t kPointAngleCut = 0.9996;
+ const Float_t kCausality0Cut = 0.19;
+ const Float_t kLikelihood01Cut = 0.45;
+ const Float_t kLikelihood1Cut = 0.5;
+ const Float_t kCombinedCut = 0.55;
//
- Int_t entries = fTrackHypothesys.GetEntriesFast();
- for (Int_t i=0;i<entries;i++){
- TObjArray * array = (TObjArray*)fTrackHypothesys.At(i);
- if (!array) continue;
- AliITStrackMI * track = (AliITStrackMI*)array->At(fBestTrackIndex[i]);
- if (track){
- dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]);
- Float_t pt1 = TMath::Abs(track->fP4*track->GetConvConst());
- mindist[i] = TMath::Min(0.02+0.03*pt1,0.15);
- if (mindist[i]<0.05) mindist[i]=0.05;
- trackarray.AddAt(track,i);
- new (&helixes[i]) AliHelix(*track);
+ //
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ Int_t ntracks = event->GetNumberOfTracks();
+ Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
+ fOriginal.Expand(ntracks);
+ fTrackHypothesys.Expand(ntracks);
+ fBestHypothesys.Expand(ntracks);
+ //
+ AliHelix * helixes = new AliHelix[ntracks+2];
+ TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain
+ TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain
+ TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain
+ Bool_t * forbidden = new Bool_t [ntracks+2];
+ Int_t *itsmap = new Int_t [ntracks+2];
+ Float_t *dist = new Float_t[ntracks+2];
+ Float_t *normdist0 = new Float_t[ntracks+2];
+ Float_t *normdist1 = new Float_t[ntracks+2];
+ Float_t *normdist = new Float_t[ntracks+2];
+ Float_t *norm = new Float_t[ntracks+2];
+ Float_t *maxr = new Float_t[ntracks+2];
+ Float_t *minr = new Float_t[ntracks+2];
+ Float_t *minPointAngle= new Float_t[ntracks+2];
+ //
+ AliESDV0MI *pvertex = new AliESDV0MI;
+ AliITStrackMI * dummy= new AliITStrackMI;
+ dummy->SetLabel(0);
+ AliITStrackMI trackat0; //temporary track for DCA calculation
+ //
+ Float_t primvertex[3]={GetX(),GetY(),GetZ()};
+ //
+ // make its - esd map
+ //
+ for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
+ itsmap[itrack] = -1;
+ forbidden[itrack] = kFALSE;
+ maxr[itrack] = kMaxR;
+ minr[itrack] = kMinR;
+ minPointAngle[itrack] = kMinPointAngle;
+ }
+ for (Int_t itrack=0;itrack<nitstracks;itrack++){
+ AliITStrackMI * original = (AliITStrackMI*)(fOriginal.At(itrack));
+ Int_t esdindex = original->fESDtrack->GetID();
+ itsmap[esdindex] = itrack;
+ }
+ //
+ // create its tracks from esd tracks if not done before
+ //
+ for (Int_t itrack=0;itrack<ntracks;itrack++){
+ if (itsmap[itrack]>=0) continue;
+ AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
+ tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
+ tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
+ if (tpctrack->fD[0]<20 && tpctrack->fD[1]<20){
+ // tracks which can reach inner part of ITS
+ // propagate track to outer its volume - with correction for material
+ CorrectForDeadZoneMaterial(tpctrack);
}
+ itsmap[itrack] = nitstracks;
+ fOriginal.AddAt(tpctrack,nitstracks);
+ nitstracks++;
}
- // Int_t multifound=0;
- Int_t vertexall =0;
- AliESDV0MI tempvertex;
- Float_t primvertex[3]={GetX(),GetY(),GetZ()};
-
- for (Int_t itrack0=0;itrack0<entries;itrack0++){
- AliITStrackMI * track0 = (AliITStrackMI*)trackarray.At(itrack0);
- if (!track0) continue;
- if (track0->fP4>0) continue;
- if (dist[itrack0]<mindist[itrack0]) continue;
- Int_t vertexes =0;
+ //
+ // fill temporary arrays
+ //
+ for (Int_t itrack=0;itrack<ntracks;itrack++){
+ AliESDtrack * esdtrack = event->GetTrack(itrack);
+ Int_t itsindex = itsmap[itrack];
+ AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
+ if (!original) continue;
+ AliITStrackMI *bestConst = 0;
+ AliITStrackMI *bestLong = 0;
+ AliITStrackMI *best = 0;
+ //
//
- for (Int_t itrack1=0;itrack1<entries;itrack1++){
- AliITStrackMI * track1 = (AliITStrackMI*)trackarray.At(itrack1);
- if (!track1) continue;
- if (track1->fP4<0) continue;
- if (dist[itrack1]<mindist[itrack1]) continue;
+ TObjArray * array = (TObjArray*) fTrackHypothesys.At(itsindex);
+ Int_t hentries = (array==0) ? 0 : array->GetEntriesFast();
+ // Get best track with vertex constrain
+ for (Int_t ih=0;ih<hentries;ih++){
+ AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
+ if (!trackh->fConstrain) continue;
+ if (!bestConst) bestConst = trackh;
+ if (trackh->fN>5.0){
+ bestConst = trackh; // full track - with minimal chi2
+ break;
+ }
+ if (trackh->fN+trackh->fNDeadZone<=bestConst->fN+bestConst->fNDeadZone) continue;
+ bestConst = trackh;
+ break;
+ }
+ // Get best long track without vertex constrain and best track without vertex constrain
+ for (Int_t ih=0;ih<hentries;ih++){
+ AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
+ if (trackh->fConstrain) continue;
+ if (!best) best = trackh;
+ if (!bestLong) bestLong = trackh;
+ if (trackh->fN>5.0){
+ bestLong = trackh; // full track - with minimal chi2
+ break;
+ }
+ if (trackh->fN+trackh->fNDeadZone<=bestLong->fN+bestLong->fNDeadZone) continue;
+ bestLong = trackh;
+ }
+ if (!best) {
+ best = original;
+ bestLong = original;
+ }
+ trackat0 = *bestLong;
+ Double_t xx,yy,zz,alpha;
+ bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
+ alpha = TMath::ATan2(yy,xx);
+ trackat0.Propagate(alpha,0);
+ // calculate normalized distances to the vertex
+ //
+ Float_t ptfac = (1.+100.*TMath::Abs(trackat0.fP4));
+ if ( bestLong->fN>3 ){
+ dist[itsindex] = trackat0.fP0;
+ norm[itsindex] = ptfac*TMath::Sqrt(trackat0.fC00);
+ normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]);
+ normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.fC11)));
+ normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+ if (!bestConst){
+ if (bestLong->fN+bestLong->fNDeadZone<6) normdist[itsindex]*=2.;
+ if (bestLong->fN+bestLong->fNDeadZone<5) normdist[itsindex]*=2.;
+ if (bestLong->fN+bestLong->fNDeadZone<4) normdist[itsindex]*=2.;
+ }else{
+ if (bestConst->fN+bestConst->fNDeadZone<6) normdist[itsindex]*=1.5;
+ if (bestConst->fN+bestConst->fNDeadZone<5) normdist[itsindex]*=1.5;
+ }
+ }
+ else{
+ if (bestConst&&bestConst->fN+bestConst->fNDeadZone>4.5){
+ dist[itsindex] = bestConst->fD[0];
+ norm[itsindex] = bestConst->fDnorm[0];
+ normdist0[itsindex] = TMath::Abs(bestConst->fD[0]/norm[itsindex]);
+ normdist1[itsindex] = TMath::Abs(bestConst->fD[0]/norm[itsindex]);
+ normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+ }else{
+ dist[itsindex] = trackat0.fP0;
+ norm[itsindex] = ptfac*TMath::Sqrt(trackat0.fC00);
+ normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]);
+ normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.fC11)));
+ normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+ if (TMath::Abs(trackat0.fP3)>1.05){
+ if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
+ if (normdist[itsindex]>3) {
+ minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
+ }
+ }
+ }
+ }
+ //
+ //-----------------------------------------------------------
+ //Forbid primary track candidates -
+ //
+ //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
+ //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
+ //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
+ //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
+ //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
+ //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
+ //-----------------------------------------------------------
+ if (bestConst){
+ if (bestLong->fN<4 && bestConst->fN+bestConst->fNDeadZone>4.5) forbidden[itsindex]=kTRUE;
+ if (normdist[itsindex]<3 && bestConst->fN+bestConst->fNDeadZone>5.5) forbidden[itsindex]=kTRUE;
+ if (normdist[itsindex]<2 && bestConst->fClIndex[0]>0 && bestConst->fClIndex[1]>0 ) forbidden[itsindex]=kTRUE;
+ if (normdist[itsindex]<1 && bestConst->fClIndex[0]>0) forbidden[itsindex]=kTRUE;
+ if (normdist[itsindex]<4 && bestConst->fNormChi2[0]<2) forbidden[itsindex]=kTRUE;
+ if (normdist[itsindex]<5 && bestConst->fNormChi2[0]<1) forbidden[itsindex]=kTRUE;
+ if (bestConst->fNormChi2[0]<2.5) {
+ minPointAngle[itsindex]= 0.9999;
+ maxr[itsindex] = 10;
+ }
+ }
+ //
+ //forbid daughter kink candidates
+ //
+ if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
+ Bool_t isElectron = kTRUE;
+ Bool_t isProton = kTRUE;
+ Double_t pid[5];
+ esdtrack->GetESDpid(pid);
+ for (Int_t i=1;i<5;i++){
+ if (pid[0]<pid[i]) isElectron= kFALSE;
+ if (pid[4]<pid[i]) isProton= kFALSE;
+ }
+ if (isElectron){
+ forbidden[itsindex]=kFALSE;
+ normdist[itsindex]*=-1;
+ }
+ if (isProton){
+ if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
+ normdist[itsindex]*=-1;
+ }
+
+ //
+ // Causality cuts in TPC volume
+ //
+ if (esdtrack->GetTPCdensity(0,10) >0.6) maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
+ if (esdtrack->GetTPCdensity(10,30)>0.6) maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
+ if (esdtrack->GetTPCdensity(20,40)>0.6) maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
+ if (esdtrack->GetTPCdensity(30,50)>0.6) maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
+ //
+ if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->fN<3) minr[itsindex]=100;
+ //
+ //
+ if (kFALSE){
+ cstream<<"Track"<<
+ "Tr0.="<<best<<
+ "Tr1.="<<((bestConst)? bestConst:dummy)<<
+ "Tr2.="<<bestLong<<
+ "Tr3.="<<&trackat0<<
+ "Esd.="<<esdtrack<<
+ "Dist="<<dist[itsindex]<<
+ "ND0="<<normdist0[itsindex]<<
+ "ND1="<<normdist1[itsindex]<<
+ "ND="<<normdist[itsindex]<<
+ "Pz="<<primvertex[2]<<
+ "Forbid="<<forbidden[itsindex]<<
+ "\n";
//
- AliHelix *h1 = &helixes[itrack0];
- AliHelix *h2 = &helixes[itrack1];
- Double_t distance = TestV0(h1,h2,pvertex);
- if (distance>0.4) continue;
+ }
+ trackarray.AddAt(best,itsindex);
+ trackarrayc.AddAt(bestConst,itsindex);
+ trackarrayl.AddAt(bestLong,itsindex);
+ new (&helixes[itsindex]) AliHelix(*best);
+ }
+ //
+ //
+ //
+ // first iterration of V0 finder
+ //
+ for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
+ Int_t itrack0 = itsmap[iesd0];
+ if (forbidden[itrack0]) continue;
+ AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
+ if (!btrack0) continue;
+ if (btrack0->fP4>0) continue;
+ AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
+ //
+ for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
+ Int_t itrack1 = itsmap[iesd1];
+ if (forbidden[itrack1]) continue;
+
+ AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
+ if (!btrack1) continue;
+ if (btrack1->fP4<0) continue;
+ Bool_t isGold = kFALSE;
+ if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
+ isGold = kTRUE;
+ }
+ AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
+ AliHelix &h1 = helixes[itrack0];
+ AliHelix &h2 = helixes[itrack1];
//
- if (distance>0.3*(dist[itrack1]+dist[itrack0])) continue;
- //if (distance>0.2*dist[itrack0]) continue;
- if (pvertex->fRr<0.3*(dist[itrack1]+dist[itrack0])) continue;
- if (pvertex->fRr>27) continue;
- pvertex->SetM(*track0);
- pvertex->SetP(*track1);
- pvertex->Update(primvertex);
+ // find linear distance
+ Double_t rmin =0;
//
//
- if (pvertex->fPointAngle<0.85) continue;
+ //
+ Double_t phase[2][2],radius[2];
+ Int_t points = h1.GetRPHIintersections(h2, phase, radius);
+ if (points==0) continue;
+ Double_t delta[2]={1000000,1000000};
+ rmin = radius[0];
+ h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
+ if (points==2){
+ if (radius[1]<rmin) rmin = radius[1];
+ h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
+ }
+ rmin = TMath::Sqrt(rmin);
+ Double_t distance = 0;
+ Double_t radiusC = 0;
+ Int_t iphase = 0;
+ if (points==1 || delta[0]<delta[1]){
+ distance = TMath::Sqrt(delta[0]);
+ radiusC = TMath::Sqrt(radius[0]);
+ }else{
+ distance = TMath::Sqrt(delta[1]);
+ radiusC = TMath::Sqrt(radius[1]);
+ iphase=1;
+ }
+ if (radiusC<TMath::Max(minr[itrack0],minr[itrack1])) continue;
+ if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1])) continue;
+ Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
+ if (distance>maxDist) continue;
+ Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
+ if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
//
//
- pvertex->fLab[0]=track0->GetLabel();
- pvertex->fLab[1]=track1->GetLabel();
- pvertex->fIndex[0] = track0->GetESDtrack()->GetID();
- pvertex->fIndex[1] = track1->GetESDtrack()->GetID();
- // calculate chi2s
+ // Double_t distance = TestV0(h1,h2,pvertex,rmin);
//
- pvertex->fChi2After = 0;
- pvertex->fChi2Before = 0;
- pvertex->fNBefore=0;
- pvertex->fNAfter=0;
+ // if (distance>maxDist) continue;
+ // if (pvertex->GetRr()<kMinR) continue;
+ // if (pvertex->GetRr()>kMaxR) continue;
+ AliITStrackMI * track0=btrack0;
+ AliITStrackMI * track1=btrack1;
+ // if (pvertex->GetRr()<3.5){
+ if (radiusC<3.5){
+ //use longest tracks inside the pipe
+ track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
+ track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
+ }
+ //
+ //
+ pvertex->SetM(*track0);
+ pvertex->SetP(*track1);
+ pvertex->Update(primvertex);
+ pvertex->SetClusters(track0->fClIndex,track1->fClIndex); // register clusters
+
+ if (pvertex->GetRr()<kMinR) continue;
+ if (pvertex->GetRr()>kMaxR) continue;
+ if (pvertex->GetPointAngle()<kMinPointAngle) continue;
+ if (pvertex->GetDist2()>maxDist) continue;
+ pvertex->SetLab(0,track0->GetLabel());
+ pvertex->SetLab(1,track1->GetLabel());
+ pvertex->SetIndex(0,track0->fESDtrack->GetID());
+ pvertex->SetIndex(1,track1->fESDtrack->GetID());
- for (Int_t i=0;i<6;i++){
- Float_t radius = fgLayers[i].GetR();
- if (pvertex->fRr>radius+0.5){
- pvertex->fNBefore+=2.;
- //
- if (track0->fClIndex[i]<=0) {
- pvertex->fChi2Before+=9;
- }else{
- Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
- track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
- pvertex->fChi2Before+=chi2;
- }
+ //
+ AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
+ AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
- if (track1->fClIndex[i]<=0) {
- pvertex->fChi2Before+=9;
- }else{
- Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
- track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
- pvertex->fChi2Before+=chi2;
+ //
+ //
+ TObjArray * array0b = (TObjArray*)fBestHypothesys.At(itrack0);
+ if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->fP3)<1.1)
+ FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
+ TObjArray * array1b = (TObjArray*)fBestHypothesys.At(itrack1);
+ if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->fP3)<1.1)
+ FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
+ //
+ AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
+ AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
+ AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
+ AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
+
+ Float_t minchi2before0=16;
+ Float_t minchi2before1=16;
+ Float_t minchi2after0 =16;
+ Float_t minchi2after1 =16;
+ Int_t maxLayer = GetNearestLayer(pvertex->GetXrp());
+
+ if (array0b) for (Int_t i=0;i<5;i++){
+ // best track after vertex
+ AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
+ if (!btrack) continue;
+ if (btrack->fN>track0l->fN) track0l = btrack;
+ // if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
+ if (btrack->fX<pvertex->GetRr()-2.) {
+ if ( (maxLayer>i+2|| (i==0)) && btrack->fN==(6-i)&&i<3){
+ Float_t sumchi2= 0;
+ Float_t sumn = 0;
+ if (maxLayer<3){ // take prim vertex as additional measurement
+ if (normdist[itrack0]>0 && htrackc0){
+ sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
+ }else{
+ sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
+ }
+ sumn += 3-maxLayer;
+ }
+ for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
+ sumn+=1.;
+ if (!btrack->fClIndex[ilayer]){
+ sumchi2+=25;
+ continue;
+ }else{
+ Int_t c=( btrack->fClIndex[ilayer] & 0x0fffffff);
+ for (Int_t itrack=0;itrack<4;itrack++){
+ if (fgLayers[ilayer].fClusterTracks[itrack][c]>=0 && fgLayers[ilayer].fClusterTracks[itrack][c]!=itrack0){
+ sumchi2+=18.; //shared cluster
+ break;
+ }
+ }
+ sumchi2+=btrack->fDy[ilayer]*btrack->fDy[ilayer]/(btrack->fSigmaY[ilayer]*btrack->fSigmaY[ilayer]);
+ sumchi2+=btrack->fDz[ilayer]*btrack->fDz[ilayer]/(btrack->fSigmaZ[ilayer]*btrack->fSigmaZ[ilayer]);
+ }
+ }
+ sumchi2/=sumn;
+ if (sumchi2<minchi2before0) minchi2before0=sumchi2;
}
+ continue; //safety space - Geo manager will give exact layer
}
-
- if (pvertex->fRr<radius-0.5){
- pvertex->fNAfter+=2.;
- //
- if (track0->fClIndex[i]<=0) {
- pvertex->fChi2After+=9;
- }else{
- Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
- track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
- pvertex->fChi2After+=chi2;
- }
-
- if (track1->fClIndex[i]<=0) {
- pvertex->fChi2After+=9;
- }else{
- Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
- track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
- pvertex->fChi2After+=chi2;
+ track0b = btrack;
+ minchi2after0 = btrack->fNormChi2[i];
+ break;
+ }
+ if (array1b) for (Int_t i=0;i<5;i++){
+ // best track after vertex
+ AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
+ if (!btrack) continue;
+ if (btrack->fN>track1l->fN) track1l = btrack;
+ // if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
+ if (btrack->fX<pvertex->GetRr()-2){
+ if ((maxLayer>i+2 || (i==0))&&btrack->fN==(6-i)&&(i<3)){
+ Float_t sumchi2= 0;
+ Float_t sumn = 0;
+ if (maxLayer<3){ // take prim vertex as additional measurement
+ if (normdist[itrack1]>0 && htrackc1){
+ sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
+ }else{
+ sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
+ }
+ sumn += 3-maxLayer;
+ }
+ for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
+ sumn+=1.;
+ if (!btrack->fClIndex[ilayer]){
+ sumchi2+=30;
+ continue;
+ }else{
+ Int_t c=( btrack->fClIndex[ilayer] & 0x0fffffff);
+ for (Int_t itrack=0;itrack<4;itrack++){
+ if (fgLayers[ilayer].fClusterTracks[itrack][c]>=0 && fgLayers[ilayer].fClusterTracks[itrack][c]!=itrack1){
+ sumchi2+=18.; //shared cluster
+ break;
+ }
+ }
+ sumchi2+=btrack->fDy[ilayer]*btrack->fDy[ilayer]/(btrack->fSigmaY[ilayer]*btrack->fSigmaY[ilayer]);
+ sumchi2+=btrack->fDz[ilayer]*btrack->fDz[ilayer]/(btrack->fSigmaZ[ilayer]*btrack->fSigmaZ[ilayer]);
+ }
+ }
+ sumchi2/=sumn;
+ if (sumchi2<minchi2before1) minchi2before1=sumchi2;
}
+ continue; //safety space - Geo manager will give exact layer
}
+ track1b = btrack;
+ minchi2after1 = btrack->fNormChi2[i];
+ break;
+ }
+ //
+ // position resolution - used for DCA cut
+ Float_t sigmad = track0b->fC00+track0b->fC11+track1b->fC00+track1b->fC11+
+ (track0b->fX-pvertex->GetRr())*(track0b->fX-pvertex->GetRr())*(track0b->fC22+track0b->fC33)+
+ (track1b->fX-pvertex->GetRr())*(track1b->fX-pvertex->GetRr())*(track1b->fC22+track1b->fC33);
+ sigmad =TMath::Sqrt(sigmad)+0.04;
+ if (pvertex->GetRr()>50){
+ Double_t cov0[15],cov1[15];
+ track0b->fESDtrack->GetInnerExternalCovariance(cov0);
+ track1b->fESDtrack->GetInnerExternalCovariance(cov1);
+ sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
+ (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
+ (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
+ sigmad =TMath::Sqrt(sigmad)+0.05;
}
- if (pvertex->fNBefore>2){
- if (pvertex->fChi2Before/pvertex->fNBefore<4.) continue; //have clusters before vetex
+ //
+ AliESDV0MI vertex2;
+ vertex2.SetM(*track0b);
+ vertex2.SetP(*track1b);
+ vertex2.Update(primvertex);
+ if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetPointAngle()>=pvertex->GetPointAngle())){
+ pvertex->SetM(*track0b);
+ pvertex->SetP(*track1b);
+ pvertex->Update(primvertex);
+ pvertex->SetClusters(track0b->fClIndex,track1b->fClIndex); // register clusters
+ pvertex->SetIndex(0,track0->fESDtrack->GetID());
+ pvertex->SetIndex(1,track1->fESDtrack->GetID());
}
- distance = FindBestPair(itrack0,itrack1,pvertex);
- if (pvertex->fPointAngle<0.85) continue;
+ pvertex->SetDistSigma(sigmad);
+ pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
+ pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
//
- if (distance>0.3) continue;
- // if (pvertex->fDistSigma*6>pvertex->fRr) continue;
- // if (pvertex->fDistSigma>0.4) continue;
- //if (pvertex->fDistNorm>5.5) continue;
- new (&oldvertexarray[vertexall]) AliV0vertex(*track0,*track1) ;
- vertexes++;
- vertexall++;
- pvertex = &vertexarray[vertexall];
- }
- }
- // printf("\n\n\nMultifound\t%d\n\n\n",multifound);
- //
- // sort vertices according quality
- Float_t quality[10000];
- Int_t indexes[10000];
- Int_t trackvertices[30000];
- for (Int_t i=0;i<entries;i++) trackvertices[i]=0;
- for (Int_t i=0;i<vertexall;i++) quality[i] = 1.-vertexarray[i].fPointAngle;
- TMath::Sort(vertexall,quality,indexes,kFALSE);
-
- for (Int_t i=0;i<vertexall;i++){
- pvertex= &vertexarray[indexes[i]];
- Int_t index0 = vertexarray[indexes[i]].fIndex[0];
- Int_t index1 = vertexarray[indexes[i]].fIndex[1];
- vertexarray[indexes[i]].fOrder[2] = i;
- vertexarray[indexes[i]].fOrder[1] = trackvertices[index1];
- vertexarray[indexes[i]].fOrder[0] = trackvertices[index0];
- if (trackvertices[index1]+trackvertices[index0]>2) continue;
- if (trackvertices[index1]+trackvertices[index0]>0) {
- // if (pvertex->fPointAngle<0.995) continue;
- }
- trackvertices[index0]++;
- trackvertices[index1]++;
-
- // event->AddV0(&oldvertexarray[indexes[i]]);
- Int_t v0index = event->AddV0MI(&vertexarray[indexes[i]]);
- AliESDtrack * ptrack0 = event->GetTrack(vertexarray[indexes[i]].fIndex[0]);
- AliESDtrack * ptrack1 = event->GetTrack(vertexarray[indexes[i]].fIndex[1]);
- if (!ptrack0 || !ptrack1){
- printf("BBBBBBBUUUUUUUUUUGGGGGGGGGG\n");
- }
- Int_t v0index0[3]={ptrack0->GetV0Index(0),ptrack0->GetV0Index(1),ptrack0->GetV0Index(2)};
- Int_t v0index1[3]={ptrack1->GetV0Index(0),ptrack1->GetV0Index(1),ptrack1->GetV0Index(2)};
- for (Int_t i=0;i<3;i++){
- if (v0index0[i]<0) {
- v0index0[i]=v0index;
- ptrack0->SetV0Indexes(v0index0);
- break;
+ // define likelihhod and causalities
+ //
+ Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
+ if (maxLayer<1){
+ Float_t fnorm0 = normdist[itrack0];
+ if (fnorm0<0) fnorm0*=-3;
+ Float_t fnorm1 = normdist[itrack1];
+ if (fnorm1<0) fnorm1*=-3;
+ if (pvertex->GetAnglep()[2]>0.1 || (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 || pvertex->GetRr()<3){
+ pb0 = TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
+ pb1 = TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
+ }
+ pvertex->SetChi2Before(normdist[itrack0]);
+ pvertex->SetChi2After(normdist[itrack1]);
+ pvertex->SetNAfter(0);
+ pvertex->SetNBefore(0);
+ }else{
+ pvertex->SetChi2Before(minchi2before0);
+ pvertex->SetChi2After(minchi2before1);
+ if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
+ pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
+ pb1 = TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
+ }
+ pvertex->SetNAfter(maxLayer);
+ pvertex->SetNBefore(maxLayer);
}
- }
- for (Int_t i=0;i<3;i++){
- if (v0index1[i]<0) {
- v0index1[i]=v0index;
- ptrack1->SetV0Indexes(v0index1);
- break;
+ if (pvertex->GetRr()<90){
+ pa0 *= TMath::Min(track0->fESDtrack->GetTPCdensity(0,60),Float_t(1.));
+ pa1 *= TMath::Min(track1->fESDtrack->GetTPCdensity(0,60),Float_t(1.));
+ }
+ if (pvertex->GetRr()<20){
+ pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
+ pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
+ }
+ //
+ pvertex->SetCausality(pb0,pb1,pa0,pa1);
+ //
+ // Likelihood calculations - apply cuts
+ //
+ Bool_t v0OK = kTRUE;
+ Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
+ p12 += pvertex->GetParamM()->GetParameter()[4]*pvertex->GetParamM()->GetParameter()[4];
+ p12 = TMath::Sqrt(p12); // "mean" momenta
+ Float_t sigmap0 = 0.0001+0.001/(0.1+pvertex->GetRr());
+ Float_t sigmap = 0.5*sigmap0*(0.6+0.4*p12); // "resolution: of point angle - as a function of radius and momenta
+
+ Float_t causalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
+ Float_t causalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Float_t(0.7))*
+ TMath::Min(pvertex->GetCausalityP()[3],Float_t(0.7)));
+ //
+ Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
+
+ Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetPointAngle())/sigmap)+
+ 0.4*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/(4.*sigmap))+
+ 0.4*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/(8.*sigmap))+
+ 0.1*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/0.01);
+ //
+ if (causalityA<kCausality0Cut) v0OK = kFALSE;
+ if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut) v0OK = kFALSE;
+ if (likelihood1<kLikelihood1Cut) v0OK = kFALSE;
+ if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
+
+ //
+ //
+ if (kFALSE){
+ Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
+ cstream<<"It0"<<
+ "Tr0.="<<track0<< //best without constrain
+ "Tr1.="<<track1<< //best without constrain
+ "Tr0B.="<<track0b<< //best without constrain after vertex
+ "Tr1B.="<<track1b<< //best without constrain after vertex
+ "Tr0C.="<<htrackc0<< //best with constrain if exist
+ "Tr1C.="<<htrackc1<< //best with constrain if exist
+ "Tr0L.="<<track0l<< //longest best
+ "Tr1L.="<<track1l<< //longest best
+ "Esd0.="<<track0->fESDtrack<< // esd track0 params
+ "Esd1.="<<track1->fESDtrack<< // esd track1 params
+ "V0.="<<pvertex<< //vertex properties
+ "V0b.="<<&vertex2<< //vertex properties at "best" track
+ "ND0="<<normdist[itrack0]<< //normalize distance for track0
+ "ND1="<<normdist[itrack1]<< //normalize distance for track1
+ "Gold="<<gold<< //
+ // "RejectBase="<<rejectBase<< //rejection in First itteration
+ "OK="<<v0OK<<
+ "rmin="<<rmin<<
+ "sigmad="<<sigmad<<
+ "\n";
+ }
+ //if (rejectBase) continue;
+ //
+ pvertex->SetStatus(0);
+ // if (rejectBase) {
+ // pvertex->SetStatus(-100);
+ //}
+ if (pvertex->GetPointAngle()>kMinPointAngle2) {
+ pvertex->SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
+ if (v0OK){
+ // AliV0vertex vertexjuri(*track0,*track1);
+ // vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
+ // event->AddV0(&vertexjuri);
+ pvertex->SetStatus(100);
+ }
+ event->AddV0MI(pvertex);
}
}
}
- delete[] vertexarray;
- delete[] oldvertexarray;
-}
-
-
-Double_t AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1, AliESDV0MI *vertex)
-{
- //
- // try to find best pair from the tree of track hyp.
- //
- TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(esdtrack0);
- Int_t entries0 = array0->GetEntriesFast();
- TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(esdtrack1);
- Int_t entries1 = array1->GetEntriesFast();
- // Float_t primvertex[3]={GetX(),GetY(),GetZ()};
- //
- //
- //AliESDV0MI v0;
- Double_t criticalradius = vertex->fRr;
//
- AliITStrackMI * track0= (AliITStrackMI*)array0->At(fBestTrackIndex[esdtrack0]);
- AliITStrackMI * track1= (AliITStrackMI*)array1->At(fBestTrackIndex[esdtrack1]);
//
- // find the best tracks after decay point
- for (Int_t itrack0=0;itrack0<entries0;itrack0++){
- AliITStrackMI * track = (AliITStrackMI*)array0->At(itrack0);
- if (!track) continue;
- if (track->fX<criticalradius-1) continue;
- if (track->fX>criticalradius){
- track0 = track;
- break;
- }
- }
-
- for (Int_t itrack1=0;itrack1<entries1;itrack1++){
- AliITStrackMI * track = (AliITStrackMI*)array1->At(itrack1);
- if (!track) continue;
- if (track->fX<criticalradius-1) continue;
- if (track->fX>criticalradius){
- track1 = track;
- break;
- }
- }
- //propagate to vertex
- Double_t alpha = TMath::ATan2(vertex->fXr[1],vertex->fXr[0]);
- Double_t radius =TMath::Sqrt(vertex->fXr[1]*vertex->fXr[1]+vertex->fXr[0]*vertex->fXr[0]);
- AliITStrackMI track0p = *track0;
- AliITStrackMI track1p = *track1;
- if (!track0p.Propagate(alpha,radius)) return 100000;
- if (!track1p.Propagate(alpha,radius)) return 100000;
- //track0p.Propagate(alpha,radius);
- //track1p.Propagate(alpha,radius);
- //
- //
- //vertex2->SetM(track0p);
- //vertex2->SetP(track1p);
- //vertex2->Update(primvertex);
- //
- AliHelix h0(track0p);
- AliHelix h1(track1p);
- Double_t distance = TestV0(&h0,&h1,vertex);
- Float_t v[3]={GetX(),GetY(),GetZ()};
- // vertex->Update(v,track0->fMass, track1->fMass);
- vertex->SetM(*track0);
- vertex->SetP(*track1);
- vertex->Update(v);
- Float_t sigma = TMath::Sqrt(track1p.GetSigmaY2()+track1p.GetSigmaZ2()+track0p.GetSigmaY2()+track0p.GetSigmaZ2());
- vertex->fDistNorm = distance/sigma;
- vertex->fDistSigma = sigma;
- return distance;
+ // delete temporary arrays
+ //
+ delete[] minPointAngle;
+ delete[] maxr;
+ delete[] minr;
+ delete[] norm;
+ delete[] normdist;
+ delete[] normdist1;
+ delete[] normdist0;
+ delete[] dist;
+ delete[] itsmap;
+ delete[] helixes;
+ delete pvertex;
}
-
-Double_t AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliESDV0MI *vertex)
+void AliITStrackerMI::RefitV02(AliESD *event)
{
//
- // test the helixes for the distnce calculate vertex characteristic
- //
- Float_t distance1,distance2;
- AliHelix & dhelix1 = *helix1;
- Double_t pp[3],xx[3];
- dhelix1.GetMomentum(0,pp,0);
- dhelix1.Evaluate(0,xx);
- AliHelix &mhelix = *helix2;
- //
- //find intersection linear
- //
- Double_t phase[2][2],radius[2];
- Int_t points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
- Double_t delta1=10000,delta2=10000;
-
- if (points>0){
- dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
- dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
- dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
- }
- if (points==2){
- dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
- dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
- dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
- }
- distance1 = TMath::Min(delta1,delta2);
- vertex->fDist1 = TMath::Sqrt(distance1);
+ //try to refit V0s in the third path of the reconstruction
//
- //find intersection parabolic
+ TTreeSRedirector &cstream = *fDebugStreamer;
//
- points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
- delta1=10000,delta2=10000;
-
- if (points>0){
- dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
- }
- if (points==2){
- dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
- }
-
- distance2 = TMath::Min(delta1,delta2);
- vertex->fDist2 = TMath::Sqrt(distance2);
- if (distance2<100){
- if (delta1<delta2){
- //get V0 info
- dhelix1.Evaluate(phase[0][0],vertex->fXr);
- dhelix1.GetMomentum(phase[0][0],vertex->fPP);
- mhelix.GetMomentum(phase[0][1],vertex->fPM);
- dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->fAngle);
- vertex->fRr = TMath::Sqrt(radius[0]);
+ Int_t nv0s = event->GetNumberOfV0MIs();
+ Float_t primvertex[3]={GetX(),GetY(),GetZ()};
+ AliESDV0MI v0temp;
+ for (Int_t iv0 = 0; iv0<nv0s;iv0++){
+ AliESDV0MI * v0mi = event->GetV0MI(iv0);
+ if (!v0mi) continue;
+ Int_t itrack0 = v0mi->GetIndex(0);
+ Int_t itrack1 = v0mi->GetIndex(1);
+ AliESDtrack *esd0 = event->GetTrack(itrack0);
+ AliESDtrack *esd1 = event->GetTrack(itrack1);
+ if (!esd0||!esd1) continue;
+ AliITStrackMI tpc0(*esd0);
+ AliITStrackMI tpc1(*esd1);
+ Double_t alpha =TMath::ATan2(v0mi->GetXr(1),v0mi->GetXr(0));
+ if (v0mi->GetRr()>85){
+ if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
+ v0temp.SetM(tpc0);
+ v0temp.SetP(tpc1);
+ v0temp.Update(primvertex);
+ if (kFALSE) cstream<<"Refit"<<
+ "V0.="<<v0mi<<
+ "V0refit.="<<&v0temp<<
+ "Tr0.="<<&tpc0<<
+ "Tr1.="<<&tpc1<<
+ "\n";
+ if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetPointAngle()>v0mi->GetPointAngle()){
+ v0mi->SetM(tpc0);
+ v0mi->SetP(tpc1);
+ v0mi->Update(primvertex);
+ }
+ }
+ continue;
}
- else{
- dhelix1.Evaluate(phase[1][0],vertex->fXr);
- dhelix1.GetMomentum(phase[1][0],vertex->fPP);
- mhelix.GetMomentum(phase[1][1],vertex->fPM);
- dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->fAngle);
- vertex->fRr = TMath::Sqrt(radius[1]);
+ if (v0mi->GetRr()>35){
+ CorrectForDeadZoneMaterial(&tpc0);
+ CorrectForDeadZoneMaterial(&tpc1);
+ if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
+ v0temp.SetM(tpc0);
+ v0temp.SetP(tpc1);
+ v0temp.Update(primvertex);
+ if (kFALSE) cstream<<"Refit"<<
+ "V0.="<<v0mi<<
+ "V0refit.="<<&v0temp<<
+ "Tr0.="<<&tpc0<<
+ "Tr1.="<<&tpc1<<
+ "\n";
+ if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetPointAngle()>v0mi->GetPointAngle()){
+ v0mi->SetM(tpc0);
+ v0mi->SetP(tpc1);
+ v0mi->Update(primvertex);
+ }
+ }
+ continue;
}
+ CorrectForDeadZoneMaterial(&tpc0);
+ CorrectForDeadZoneMaterial(&tpc1);
+ // if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
+ if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
+ v0temp.SetM(tpc0);
+ v0temp.SetP(tpc1);
+ v0temp.Update(primvertex);
+ if (kFALSE) cstream<<"Refit"<<
+ "V0.="<<v0mi<<
+ "V0refit.="<<&v0temp<<
+ "Tr0.="<<&tpc0<<
+ "Tr1.="<<&tpc1<<
+ "\n";
+ if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetPointAngle()>v0mi->GetPointAngle()){
+ v0mi->SetM(tpc0);
+ v0mi->SetP(tpc1);
+ v0mi->Update(primvertex);
+ }
+ }
}
- //
- //
- return vertex->fDist2;
}