}
fESDtrack=&t;
+ // if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
- if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
}
void AliITStrackV2::UpdateESDtrack(ULong_t flags) const {
//------------------------------------------------------------------
Double_t x1=fX, x2=xk, dx=x2-x1;
Double_t f1=fP2, f2=f1 + fP4*dx;
- if (TMath::Abs(f2) >= 0.9999) {
- Int_t n=GetNumberOfClusters();
- if (n>kWARN)
- Warning("PropagateTo","Propagation failed !\n",n);
+ if (TMath::Abs(f2) >= 0.98) {
+ // MI change - don't propagate highly inclined tracks
+ // covariance matrix distorted
+ // Int_t n=GetNumberOfClusters();
+ // if (n>kWARN)
+ // Warning("PropagateTo","Propagation failed !\n",n);
return 0;
}
{
Double_t dx=xk-fX;
Double_t f1=fP2, f2=f1 + fP4*dx;
- if (TMath::Abs(f2) >= 0.9999) {
- Int_t n=GetNumberOfClusters();
- if (n>kWARN)
- Warning("Propagate","Propagation failed (%d) !\n",n);
+ if (TMath::Abs(f2) >= 0.98) {
+ // don't propagate highly inclined tracks MI
return 0;
}
+ // if (TMath::Abs(f2) >= 0.9999) {
+// Int_t n=GetNumberOfClusters();
+// if (n>kWARN)
+// Warning("Propagate","Propagation failed (%d) !\n",n);
+// return 0;
+// }
Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
fX=xk;
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) {
//--------------------------------------------------------------------
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);
// 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.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
//
}
//GetBestHypothesysMIP(itsTracks);
- FindV0(event);
+ UpdateTPCV0(event);
+ FindV02(event);
fAfterV0 = kTRUE;
//GetBestHypothesysMIP(itsTracks);
//
}
fTrackHypothesys.Delete();
+ fBestHypothesys.Delete();
fOriginal.Clear();
delete []fCoeficients;
fCoeficients=0;
if (fTrackToFollow.Propagate(fv+a,xv)) {
fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
- Float_t d=fTrackToFollow.GetD(GetX(),GetY());
- Float_t z=fTrackToFollow.GetZ()-GetZ();
- fTrackToFollow.GetESDtrack()->SetImpactParameters(d,z);
//UseClusters(&fTrackToFollow);
{
AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ());
}
-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 (fConstraint[fPass]&&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;
ntracks[ilayer]++;
}
- if (fConstraint[fPass]&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){ //big theta -- for low mult. runs
+ 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;
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 = fConstraint[fPass]? 20: 5;
+ Int_t max = constrain? 20: 5;
for (Int_t i=0;i<TMath::Min(max,ntracks[0]);i++) {
AliITStrackMI & track= tracks[0][nindexes[0][i]];
if (track.GetNumberOfClusters()<2) continue;
- if (!fConstraint[fPass]&&track.fNormChi2[0]>7.)continue;
+ if (!constrain&&track.fNormChi2[0]>7.)continue;
AddTrackHypothesys(new AliITStrackMI(track), esdindex);
}
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]){
+ if (!constrain){
for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
AliITStrackMI & track= tracks[2][nindexes[2][i]];
if (track.GetNumberOfClusters()<3) continue;
- if (!fConstraint[fPass]&&track.fNormChi2[2]>7.)continue;
- if (fConstraint[fPass]) track.fNSkipped+=2;
- if (!fConstraint[fPass]){
+ 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 (!backtrack->Improve(0,xyzv,ersv)) continue;
if (!backtrack->PropagateToVertex()) continue;
backtrack->ResetCovariance();
- if (!backtrack->Improve(0,xyzv,ersv)) continue;
+ if (!backtrack->Improve(0,xyzv,ersv)) continue;
}else{
backtrack->ResetCovariance();
}
//
sigmarfi = 0.004+1.4 *TMath::Abs(track->fP4)+332.*track->fP4*track->fP4;
sigmaz = 0.011+4.37*TMath::Abs(track->fP4);
-
}
+Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
+ //
+ //Get nearest upper layer close to the point xr.
+ // rough approximation
+ //
+ const Float_t radiuses[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<radiuses[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
+ //
+ // 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;
+
+ //
+ //
+ 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
//
- //TTreeSRedirector cstream("itsv0.root");
- Int_t centries=0;
- AliHelix * helixes = new AliHelix[30000];
- TObjArray trackarray(30000);
- TObjArray trackarrayc(30000);
- Float_t * dist = new Float_t[30000];
- Float_t * normdist0 = new Float_t[30000];
- Float_t * normdist1 = new Float_t[30000];
- Float_t * normdist = new Float_t[30000];
- Float_t * norm = new Float_t[30000];
- AliESDV0MI *vertexarray = new AliESDV0MI[100000];
- AliESDV0MI *pvertex = &vertexarray[0];
- AliITStrackMI * dummy=0;
+ Float_t primvertex[3]={GetX(),GetY(),GetZ()};
//
+ // make its - esd map
//
- Int_t entries = fTrackHypothesys.GetEntriesFast();
- for (Int_t i=0;i<entries;i++){
- TObjArray * array = (TObjArray*)fTrackHypothesys.At(i);
- if (!array) continue;
- // get best track without vertex constrain
- Int_t hentries = array->GetEntriesFast();
+ 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++;
+ }
+ //
+ // 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;
+ //
//
- // best with vertex constrain
- AliITStrackMI * trackc = (AliITStrackMI*)array->At(0);
- if (trackc&&trackc->fConstrain&&trackc->fN==6&&trackc->fNormChi2[0]<2.) continue;
- trackc=0;
+ 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 (trackh->fN<6) continue;
- trackc = trackh;
- if (!dummy) dummy = trackc;
- dummy->SetLabel(0);
+ 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;
- }
- //
- // best without vertex
- AliITStrackMI * track = 0;
+ }
+ // 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;
- track = trackh;
- break;
+ 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 (trackc&&track){
- if (trackc->fChi2MIP[1]<2.) continue;
- if (trackc->fChi2MIP[0]<2. && trackc->fChi2MIP[1]<2.) continue;
- trackarrayc.AddAt(trackc,i);
- if (trackc->fN==6&&track->fN&&trackc->fNormChi2[0] < track->fNormChi2[0]-2) continue;
+ 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
//
+ if ( bestLong->fN>3 ){
+ dist[itsindex] = trackat0.fP0;
+ norm[itsindex] = TMath::Sqrt(trackat0.fC00);
+ normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]);
+ normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/TMath::Sqrt(trackat0.fC11));
+ normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+ }
+ 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] = TMath::Sqrt(trackat0.fC00);
+ normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]);
+ normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/TMath::Sqrt(trackat0.fC11));
+ normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+ if (TMath::Abs(trackat0.fP3)>1.1){
+ if (normdist[itsindex]<10) forbidden[itsindex]=kTRUE;
+ }
+ }
+ }
//
+ //-----------------------------------------------------------
+ //Forbid primary track candidates -
//
- if (track){
- dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]);
- norm[i] = track->fDnorm[0];
- normdist0[i] = TMath::Abs(track->fD[0]/track->fDnorm[0]);
- normdist1[i] = TMath::Abs(track->fD[1]/track->fDnorm[1]);
- normdist[i] = TMath::Sqrt(normdist0[i]*normdist0[i]+normdist1[i]*normdist1[i]);
- if (track->IsGoldPrimary()) continue; //primary track
- if (track->fD[0]<0.02 && (track->fN+track->fNDeadZone>5.8)){
- if (normdist[i]<3.) continue; // primary track - cutoff 3 sigma
- if (normdist0[i]<2.) continue; //DCA normalized cut 2 sigma
+ //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;
}
- trackarray.AddAt(track,i);
- new (&helixes[i]) AliHelix(*track);
}
+ //
+ //forbid daughter kink candidates
+ //
+ if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
+ Bool_t isElectron = kTRUE;
+ Double_t pid[5];
+ esdtrack->GetESDpid(pid);
+ for (Int_t i=1;i<5;i++){
+ if (pid[0]<pid[i]) isElectron= kFALSE;
+ }
+ if (isElectron){
+ forbidden[itsindex]=kFALSE;
+ normdist[itsindex]+=4;
+ }
+ //
+ // 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 (0){
+ 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";
+ //
+ }
+ trackarray.AddAt(best,itsindex);
+ trackarrayc.AddAt(bestConst,itsindex);
+ trackarrayl.AddAt(bestLong,itsindex);
+ new (&helixes[itsindex]) AliHelix(*best);
}
//
//
- // 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;
+ //
+ // 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);
//
- TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(itrack0);
- //
- Int_t vertexes =0;
- for (Int_t itrack1=0;itrack1<entries;itrack1++){
- AliITStrackMI * track1 = (AliITStrackMI*)trackarray.At(itrack1);
- if (!track1) continue;
- if (track1->fP4<0) continue;
- AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
- if (trackc0&&trackc1){
- if (TMath::Min(trackc0->fChi2MIP[1],trackc1->fChi2MIP[1])<2.) continue;
+ 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;
}
- if (track1->fNDeadZone+track0->fNDeadZone>1.1) continue;
- TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(itrack1);
+ AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
+ AliHelix &h1 = helixes[itrack0];
+ AliHelix &h2 = helixes[itrack1];
//
- //if (normdist0[itrack0]+normdist0[itrack1]<3) continue;
- //if (normdist[itrack0]+normdist[itrack1]<4) continue;
+ // find linear distance
+ Double_t rmin =0;
//
//
- AliHelix *h1 = &helixes[itrack0];
- AliHelix *h2 = &helixes[itrack1];
- Double_t rmin =0;
- Double_t distance = TestV0(h1,h2,pvertex,rmin);
//
- if (distance>0.4) continue;
- if (pvertex->GetRr()<0.3) continue;
- if (pvertex->GetRr()>20.) continue;
- pvertex->SetM(*track0);
- pvertex->SetP(*track1);
- pvertex->Update(primvertex);
- if (pvertex->GetRr()<0.3) continue;
- if (pvertex->GetRr()>20.) continue;
- if (track1->fNDeadZone+track0->fNDeadZone>0.5 &&distance>0.12) continue;
-
+ Double_t phase[2][2],radius[2];
+ Int_t points = h1.GetRPHIintersections(h2, phase, radius);
+ if (points==0) continue;
+ Double_t delta[2]={1000,1000};
+ 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 (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;
//
-
- if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2
- ||(centries<5000&&(pvertex->GetPointAngle()>0.95))){
- //cstream<<"Iter0"<<track0<<track1<<pvertex<<normdist[itrack0]<<normdist[itrack1]<<"\n";
- centries++;
- }
//
+ // Double_t distance = TestV0(h1,h2,pvertex,rmin);
//
- if (pvertex->GetPointAngle()<0.85) continue;
- // if (normdist[itrack0]+normdist[itrack1]<6&&pvertex->GetPointAngle()<0.99) continue;
+ // 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);
+ 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->GetESDtrack()->GetID());
- pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
- // calculate chi2s
- //
- pvertex->SetChi2After(0);
- pvertex->SetChi2Before(0);
- pvertex->SetNBefore(0);
- pvertex->SetNAfter(0);
- for (Int_t i=0;i<6;i++){
- Float_t radius = fgLayers[i].GetR();
- if (pvertex->GetRr()>radius+0.5){
- pvertex->SetNBefore(pvertex->GetNBefore()+2.);
- //
- if (track0->fClIndex[i]<=0) {
- pvertex->SetChi2Before(pvertex->GetChi2Before()+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->SetChi2Before(pvertex->GetChi2Before()+chi2);
- }
+ //
+ AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
+ AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
- if (track1->fClIndex[i]<=0) {
- pvertex->SetChi2Before(pvertex->GetChi2Before()+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;
- pvertex->SetChi2Before(pvertex->GetChi2Before()+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) {
+ if (maxLayer>i+2 && btrack->fN==(6-i)&&i<2){
+ Float_t sumchi2= 0;
+ Float_t sumn = 0;
+ for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
+ sumn+=1.;
+ if (!btrack->fClIndex[ilayer]){
+ sumchi2+=25;
+ continue;
+ }else{
+ 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->GetRr()<radius-0.5){
- pvertex->SetNAfter(pvertex->GetNAfter()+2.);
- //
- if (track0->fClIndex[i]<=0) {
- pvertex->SetChi2After(pvertex->GetChi2After()+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->SetChi2After(pvertex->GetChi2After()+chi2);
- }
-
- if (track1->fClIndex[i]<=0) {
- pvertex->SetChi2After(pvertex->GetChi2After()+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->SetChi2After(pvertex->GetChi2After()+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){
+ if (maxLayer>i+2&&btrack->fN==(6-i)&&(i<2)){
+ Float_t sumchi2= 0;
+ Float_t sumn = 0;
+ for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
+ sumn+=1.;
+ if (!btrack->fClIndex[ilayer]){
+ sumchi2+=30;
+ continue;
+ }else{
+ 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;
}
- if (pvertex->GetNBefore()>2){
- if (pvertex->GetChi2Before()/pvertex->GetNBefore()<4.) continue; //have clusters before vetex
- }
- Int_t ibest0=0,ibest1=0;
- AliITStrackMI * ntrack0 = track0;
- AliITStrackMI * ntrack1 = track1;
- //
- //
- //PH Float_t oldistance = pvertex->GetDist2();
- Bool_t improve = FindBestPair(itrack0,itrack1,pvertex,ibest0,ibest1); // try to improve vertex
- if (pvertex->GetPointAngle()<0.5) continue;
- distance = pvertex->GetDist2();
- if (improve){
- ntrack0 = (AliITStrackMI*)array0->At(ibest0);
- ntrack1 = (AliITStrackMI*)array1->At(ibest1);
- }
- Bool_t accept0 = kFALSE; // accept ==> because of pointing angle
- if (pvertex->GetPointAngle()>0.999){
- if (pvertex->GetRr()<3.5 && (ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone)<11.5) continue;
- if (pvertex->GetRr()>3.5&& pvertex->GetDistNorm()<12) accept0 = kTRUE;
- if (pvertex->GetRr()>1 && pvertex->GetDist2()<0.1 && pvertex->GetDistNorm()<12) accept0 = kTRUE;
- if (pvertex->GetPointAngle()>0.9995&&pvertex->GetRr()>5.) accept0 = kTRUE;
- }
- Bool_t reject1= kFALSE; // reject ==> bad kinematic
- //
- reject1 |= TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>1.02 ||
- TMath::Abs(ntrack0->fN-ntrack1->fN)>1.02; // cut1
- reject1 |= ntrack0->fNUsed+ntrack1->fNUsed>1.01; // cut2
- reject1 |= pvertex->GetDistNorm()>12; // cut3
- reject1 |= pvertex->GetDist2()>0.1 && improve; // cut4
- reject1 |= (TMath::Abs(ntrack0->fD[0])+TMath::Abs(ntrack1->fD[0]))/pvertex->GetDist2()<5; //cut5
- reject1 |= TMath::Abs(ntrack0->fD[0]/pvertex->GetDist2())<2 || TMath::Abs(ntrack1->fD[0]/pvertex->GetDist2())<2; //cut 6
//
- // small radii cuts
- Bool_t reject2 = kFALSE;
- if (pvertex->GetRr()<3.6){
- reject2 |= (TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>0.01); // cut7
- reject2 |= ntrack0->fNUsed+ntrack1->fNUsed>0.01; // cut8
- reject2 |= ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone<11.5; // cut9
- reject2 |= (ntrack0->fN+ntrack1->fN<11.5)&&pvertex->GetRr()<2; // cut10
- reject2 |= pvertex->GetDist2()>0.1; // cut11
- }
- //PH AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
- //PH AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
+ // 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;
+ }
+ //
+ 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->SetDistSigma(sigmad);
+ pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
//
+ // define likelihhod and causalities
//
- //
- if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2
- ||(centries<500000)){
- /*
- cstream<<"It1"<<"Tr0.="<<ntrack0<<"TR1.="<<ntrack1<<"V0.="<<pvertex<<"ND0.="<<normdist[itrack0]<<"ND1.="<<
- normdist[itrack1]<<"D.="<<distance<<"DistOld="<<oldistance<<"Imp.="<<improve<<
- "A0="<<accept0<<"R1="<<reject1<<"R2="<<reject2<<"Rmin.="<<rmin<<
- "TrC0.="<<htrackc0<<"TRC1.="<<htrackc1<<"\n";
- */
- centries++;
+ Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
+ if (maxLayer<2){
+ if (pvertex->GetAnglep()[2]>0.2){
+ pb0 = TMath::Exp(-TMath::Min(normdist[itrack0],Float_t(16.))/12.);
+ pb1 = TMath::Exp(-TMath::Min(normdist[itrack1],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.2){
+ 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);
}
-
- if (!accept0 && (reject1 || reject2)) continue;
-
-// if (distance>0.5) continue;
-// distance = pvertex->GetDist2();
-// if (pvertex->GetRr()>25 || pvertex->GetRr()<0.2) continue;
-// if (pvertex->GetRr()/pvertex->fDistSigma<1) continue;
-// if (pvertex->GetDistNorm()>10) continue;
-// if (pvertex->GetPointAngle()<0.85) continue;
-// if ((normdist[itrack0]<3||normdist[itrack1]<3)){
-// if (pvertex->GetPointAngle()<0.99||pvertex->GetDist2()>0.15) continue;
-// }
-// if (distance>0.05*(0.8+0.2*(0.5+pvertex->GetRr()))) continue;
-// if (pvertex->GetRr()<0.3) continue;
-// if (pvertex->GetRr()>27.) continue;
-
-
-
- //
- if (distance<0.3 &&pvertex->GetPointAngle()>0.998){
- track0->fGoldV0 = kTRUE;
- track1->fGoldV0 = kTRUE;
+ 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.));
}
- vertexes++;
- vertexall++;
- if (vertexall>=100000) break;
- pvertex = &vertexarray[vertexall];
- }
- }
- // printf("\n\n\nMultifound\t%d\n\n\n",multifound);
- //
- // sort vertices according quality
- Float_t quality[40000];
- Int_t indexes[40000];
- Int_t trackvertices[40000];
- for (Int_t i=0;i<entries;i++) trackvertices[i]=0;
- for (Int_t i=0;i<vertexall;i++) {
- Float_t norm = 1.-0.999*TMath::Abs(vertexarray[i].GetPointAngle());
- Float_t fnormdist = 1./(1+vertexarray[i].GetRr());
- quality[i] = norm*fnormdist;
- }
- //
- TMath::Sort(vertexall,quality,indexes,kFALSE);
-
- for (Int_t i=0;i<vertexall;i++){
- pvertex= &vertexarray[indexes[i]];
- Int_t index0 = vertexarray[indexes[i]].GetIndex(0);
- Int_t index1 = vertexarray[indexes[i]].GetIndex(1);
- vertexarray[indexes[i]].SetOrder(2,i);
- vertexarray[indexes[i]].SetOrder(1,trackvertices[index1]);
- vertexarray[indexes[i]].SetOrder(0,trackvertices[index0]);
- Int_t v0index = event->AddV0MI(&vertexarray[indexes[i]]);
- //
- if (trackvertices[index1]+trackvertices[index0]>5) continue;
- if (trackvertices[index0]>2) continue;
- if (trackvertices[index1]>2) continue;
-
- if (trackvertices[index1]+trackvertices[index0]>0) {
- // if (pvertex->GetPointAngle()<0.995) continue;
- }
- trackvertices[index0]++;
- trackvertices[index1]++;
-
- AliESDtrack * ptrack0 = event->GetTrack(vertexarray[indexes[i]].GetIndex(0));
- AliESDtrack * ptrack1 = event->GetTrack(vertexarray[indexes[i]].GetIndex(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;
- }
- }
- for (Int_t i=0;i<3;i++){
- if (v0index1[i]<0) {
- v0index1[i]=v0index;
- ptrack1->SetV0Indexes(v0index1);
- break;
+ 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;
}
- }
- }
-
- delete [] helixes;
- delete [] dist;
- delete [] normdist0;
- delete [] normdist1;
- delete [] normdist;
- delete [] norm;
-
- delete[] vertexarray;
-}
-
-
-
-Bool_t AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1, AliESDV0MI *vertex, Int_t &i0, Int_t &i1)
-{
- //
- // 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();
- // AliITStrackMI *orig0 = (AliITStrackMI*)fOriginal.At(esdtrack0);
- //AliITStrackMI *orig1 = (AliITStrackMI*)fOriginal.At(esdtrack1);
- Double_t criticalradius = vertex->GetRr();
- AliITStrackMI * track0= 0;
- AliITStrackMI * track1= 0;
- i0 = -1;
- i1 = -1;
- //
- //
- Float_t rfirst0[2000]; //radius position of the first cluster - track0
- Float_t rfirst1[2000]; // - track1
- Float_t maxlocalx0=0; //local x for first track
- Float_t maxlocalx1=0; //local x for second track
- Float_t cs0=1, sn0=0; //rotations
- Float_t cs1=1, sn1=0; //rotations
-
- //
- for (Int_t itrack0=0;itrack0<entries0;itrack0++){
- rfirst0[itrack0]=-1.;
- AliITStrackMI * htrack0 = (AliITStrackMI*)array0->At(itrack0);
- if (!htrack0) continue;
- if (htrack0->fConstrain) continue;
- if (i0<0){
- i0 = itrack0;
- track0 = htrack0;
- }
- Double_t cs = TMath::Cos(htrack0->fAlpha);
- Double_t sn = TMath::Sin(htrack0->fAlpha);
- Double_t x = htrack0->fX*cs - htrack0->fP0*sn;
- Double_t y = htrack0->fX*sn + htrack0->fP0*cs;
- Double_t radius = TMath::Sqrt(x*x+y*y);
- if (criticalradius<3 && radius>6&&htrack0->fNDeadZone<0.2) continue; // all cluster required
- if (criticalradius>10 && radius<6) continue; // causality
- Double_t localx = TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn);
- if (localx>maxlocalx0) {
- maxlocalx0=localx;
- cs0 = cs; sn0=sn;
- }
- rfirst0[itrack0] = radius;
- }
- for (Int_t itrack1=0;itrack1<entries1;itrack1++){
- rfirst1[itrack1]=-1.;
- AliITStrackMI * htrack1 = (AliITStrackMI*)array1->At(itrack1);
- if (!htrack1) continue;
- if (htrack1->fConstrain) continue;
- if (i1<0){
- i1 = itrack1;
- track1 = htrack1;
- }
- Double_t cs = TMath::Cos(htrack1->fAlpha);
- Double_t sn = TMath::Sin(htrack1->fAlpha);
- //
- //
- Double_t x = htrack1->fX*cs - htrack1->fP0*sn;
- Double_t y = htrack1->fX*sn + htrack1->fP0*cs;
- Double_t radius = TMath::Sqrt(x*x+y*y);
- if (criticalradius<3 && radius>6&&htrack1->fNDeadZone<0.2) continue; //all clusters required
- if (criticalradius>10 && radius<6) continue; //causality
- Double_t localx = TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn);
- if (localx>maxlocalx1) {
- maxlocalx1=localx;
- cs1 = cs; sn1=sn;
- }
- rfirst1[itrack1] = radius;
- }
- //
- //
- //
- const Float_t radiuses[4]={4,6.5,15.03,24.};
- //
- //
- // find the best tracks after decay point
- Float_t bestquality =100000;
- Float_t bestradius=0;
- AliESDV0MI v0;
- Double_t vpos[3];
- Float_t v[3]={GetX(),GetY(),GetZ()};
- //
- //
- for (Int_t itrack0=0;itrack0<entries0;itrack0++){
- if (rfirst0[itrack0]<0) continue;
- AliITStrackMI * htrack0 = (AliITStrackMI*)array0->At(itrack0);
- if (!htrack0) continue;
- //
- for (Int_t itrack1=0;itrack1<entries1;itrack1++){
- if (rfirst1[itrack1]<0) continue;
- AliITStrackMI * htrack1 = (AliITStrackMI*)array1->At(itrack1);
- if (!htrack1) continue;
- if (TMath::Abs(rfirst0[itrack0]-rfirst1[itrack1])>1.) continue;
- if (htrack0->fClIndex[6-htrack0->fN]==htrack1->fClIndex[6-htrack1->fN]) continue; //sharing of last cluster not allowe
//
- if (htrack0->fNUsed+htrack1->fNUsed>0) continue; //sharing of clusters not allowed
+ pvertex->SetCausality(pb0,pb1,pa0,pa1);
//
- //
- v0.SetM(*htrack0);
- v0.SetP(*htrack1);
- // if (v0.Update(v)==0) continue;
- v0.Update(v);
- if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) <v0.GetRr()-0.3) continue;
+ // 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 (v0.GetDist2()>0.3) continue;
- if (v0.GetRr()<radiuses[1] && ( TMath::Abs(htrack0->fN+htrack0->fNDeadZone-htrack1->fN-htrack1->fNDeadZone)>0.5)) continue;
+ 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 (v0.GetRr()<3. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<11.7) continue;
- //if (v0.GetRr()<6. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<9.7) continue;
- if (v0.GetRr()<3. && (htrack0->fN+htrack1->fN)<11.7) continue;
- if (v0.GetRr()<6. && (htrack0->fN+htrack1->fN)<9.7) continue;
- Double_t localx0=v0.GetXr(0)*cs0+v0.GetXr(1)*sn0;
- Double_t localx1=v0.GetXr(0)*cs1+v0.GetXr(1)*sn1;
- Double_t maxlocalx = TMath::Max(localx0,localx1);
- if (maxlocalx<3.4 && (htrack0->fN+htrack1->fN)<11.7) continue;
- if (maxlocalx<6.1 && (htrack0->fN+htrack1->fN)<9.7) continue;
//
- Float_t fnormdist = v0.GetDist2()/0.05;
- fnormdist +=(htrack0->fNormChi2[6-htrack0->fN]+htrack1->fNormChi2[6-htrack1->fN]);
- fnormdist +=3*(htrack0->fNUsed+htrack1->fNUsed);
- if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) <v0.GetRr()+0.2){
- fnormdist +=(v0.GetRr()+0.2-TMath::Min(rfirst0[itrack0],rfirst1[itrack1]))/0.1;
- }
+ if (kTRUE){
+ 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;
//
- Float_t quality = fnormdist;
- if (quality<bestquality && v0.GetDist2()<vertex->GetDist2()){
- i0=itrack0;
- i1=itrack1;
- track0 =htrack0;
- track1 =htrack1;
- bestquality = quality;
- bestradius = v0.GetRr();
- vpos[0] = v0.GetXr(0);
- vpos[1] = v0.GetXr(1);
- vpos[2] = v0.GetXr(2);
+ pvertex->SetStatus(0);
+ // if (rejectBase) {
+ // pvertex->SetStatus(-100);
+ //}
+ if (pvertex->GetPointAngle()>kMinPointAngle2) {
+ if (v0OK){
+ AliV0vertex vertexjuri(*track0,*track1);
+ vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
+ pvertex->SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
+ event->AddV0(&vertexjuri);
+ pvertex->SetStatus(100);
+ }
+ event->AddV0MI(pvertex);
}
}
}
//
//
- if (!track0||!track1) return kFALSE;
- if (track0->fNUsed+track1->fNUsed>0) return kFALSE; // sharing of clusters not allowed
- //
- //propagate to vertex
-
- Double_t alpha = TMath::ATan2(vpos[1],vpos[0]);
- //
- AliITStrackMI track0p = *track0;
- AliITStrackMI track1p = *track1;
- //
- // RefitAt(bestradius+0.5,&track0p,track0);
- //RefitAt(bestradius+0.5,&track1p,track1);
-
- if ( track0p.Propagate(alpha,bestradius+0.2)) {
- track0= &track0p;
- }
- if (track1p.Propagate(alpha,bestradius+0.2)){
- track1 = &track1p;
- }
- //
- v0.SetM(*track0);
- v0.SetP(*track1);
- v0.Update(v);
- if (v0.GetDist2()<vertex->GetDist2() && v0.GetRr()<20){
- vertex->SetM(*track0);
- vertex->SetP(*track1);
- vertex->Update(v);
- return kTRUE;
- }
- return kFALSE;
+ // 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, Double_t &rmin)
-{
- //
- // test the helixes for the distnce calculate vertex characteristic
- //
- rmin =0;
- 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){
- rmin = radius[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){
- if (radius[1]<rmin) rmin = radius[1];
- 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);
- }
- rmin = TMath::Sqrt(rmin);
- distance1 = TMath::Min(delta1,delta2);
- vertex->SetDist1(TMath::Sqrt(distance1));
-
- //
- //find intersection parabolic
- //
- 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->SetDist2(TMath::Sqrt(distance2));
- if (distance2<100){
- if (delta1<delta2){
- //get V0 info
- dhelix1.Evaluate(phase[0][0],vertex->GetXrp());
- dhelix1.GetMomentum(phase[0][0],vertex->GetPPp());
- mhelix.GetMomentum(phase[0][1],vertex->GetPMp());
- dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->GetAnglep());
- vertex->SetRr(TMath::Sqrt(radius[0]));
- }
- else{
- dhelix1.Evaluate(phase[1][0],vertex->GetXrp());
- dhelix1.GetMomentum(phase[1][0],vertex->GetPPp());
- mhelix.GetMomentum(phase[1][1],vertex->GetPMp());
- dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->GetAnglep());
- vertex->SetRr(TMath::Sqrt(radius[1]));
- }
- }
- //
- //
- return vertex->GetDist2();
-}
+
class AliV0vertex;
class AliESDV0MI;
+class TTreeSRedirector;
+
Int_t GetSkip() const {return fSkip;}
void SetSkip(Int_t skip){fSkip=skip;}
void IncAccepted(){fAccepted++;}
- Int_t GetAccepted() const {return fAccepted;}
+ Int_t GetAccepted() const {return fAccepted;}
protected:
AliITSlayer(const AliITSlayer& /*layer*/){;}
Double_t fR; // mean radius of this layer
AliITStrackerMI::AliITSdetector & GetDetector(Int_t layer, Int_t n) const {return GetLayer(layer).GetDetector(n); }
protected:
- void FindV0(AliESD *event); //try to find V0
- Double_t TestV0(AliHelix *h1, AliHelix *h2, AliESDV0MI *vertex, Double_t &rmin); //try to find V0 - return DCA
- Bool_t FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliESDV0MI *vertex, Int_t &ibest0, Int_t &ibest1); // try to find best pair from the tree of track hyp.
+ Int_t GetNearestLayer(const Double_t *xr) const; //get nearest upper layer close to the point xr
+ void FindV02(AliESD *event); //try to find V0
+ void UpdateTPCV0(AliESD *event); //try to update, or reject TPC V0s
void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
void CookLabel(AliITStrackMI *t,Float_t wrong) const;
Double_t GetEffectiveThickness(Double_t y, Double_t z) const;
- void FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex);
+ void FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain);
void ResetBestTrack() {
fBestTrack.~AliITStrackMI();
new(&fBestTrack) AliITStrackMI(fTrackToFollow);
AliITStrackMI fBestTrack; // "best" track
AliITStrackMI fTrackToFollow; // followed track
TObjArray fTrackHypothesys; // ! array with track hypothesys- ARRAY is the owner of tracks- MI
+ TObjArray fBestHypothesys; // ! array with track hypothesys- ARRAY is the owner of tracks- MI
TObjArray fOriginal; // ! array with seeds from the TPC
Int_t fBestTrackIndex[100000]; // ! index of the best track
Int_t fCurrentEsdTrack; // ! current esd track - MI
Int_t fLayersNotToSkip[kMaxLayer]; // layer masks
Int_t fLastLayerToTrackTo; // the innermost layer to track to
Float_t * fCoeficients; //! working array with errors and mean cluser shape
+ AliESD * fEsd; //! pointer to the ESD event
+ TTreeSRedirector *fDebugStreamer; //!debug streamer
private:
AliITStrackerMI(const AliITStrackerMI * /*tracker*/){;}
ClassDef(AliITStrackerMI,2) //ITS tracker MI
//--------------------------------------------------------------------
Float_t circ=TMath::TwoPi()*fR;
Int_t sec=Int_t(kNsector*c->GetPhiR()/circ);
+ if (sec>=kNsector) {
+ ::Error("InsertCluster","Wrong sector !\n");
+ return 1;
+ }
Int_t &n=fN[sec];
if (n>=kMaxClusterPerSector) {
::Error("InsertCluster","Too many clusters !\n");
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
+TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
+TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
+TCut crec("crec","fReconstructed==1");
+TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
+TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
+
+TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
+TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
+TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
+TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.");
+
+
+
+void MakeAliases(AliESDComparisonDraw&comp)
+{
+ //
+ // aliases definition
+ //
+ comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+ comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
+ comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
+ comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
+ comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
+ comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+ comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+
+ comp.fTree->SetAlias("trddedx","(RC.fESDTrack.fTRDsignals[0]+RC.fESDTrack.fTRDsignals[1]+RC.fESDTrack.fTRDsignals[2]+RC.fESDTrack.fTRDsignals[3]+RC.fESDTrack.fTRDsignals[4]+RC.fESDTrack.fTRDsignals[5])/6.");
+
+ comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
+ comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)");
+
+ comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]");
+ comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum");
+ comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum");
+ comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum");
+ comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum");
+ comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum");
+ comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
+}
+
+
+
void AliESDRecInfo::UpdatePoints(AliESDtrack*track)
{
}
- if (fStatus[1]>0 &&info->fNTPCRef>0){
+ if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
//TPC
fESDTrack.GetInnerXYZ(fTPCinR1);
fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
fNextTreeGenEntryToRead = 0;
fNextKinkToRead = 0;
+ fNextV0ToRead =0;
cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
eventNr++) {
fRecInfo->fESDTrack =*track;
if (track->GetITStrack())
fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
+ else{
+ fRecInfo->fITStrack = *track;
+ }
if (track->GetTRDtrack()){
fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
}
//
//Dafault constructor
//
+ for (Int_t i=0;i<4;i++){fCausality[i]=0;}
+}
+
+
+void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
+{
+ //
+ // set probabilities
+ //
+ fCausality[0] = pb0; // probability - track 0 exist before vertex
+ fCausality[1] = pb1; // probability - track 1 exist before vertex
+ fCausality[2] = pa0; // probability - track 0 exist close after vertex
+ fCausality[3] = pa1; // probability - track 1 exist close after vertex
}
void AliESDV0MI::SetP(const AliExternalTrackParam & paramp) {
//
- // set mother
+ // set track +
//
fParamP = paramp;
}
void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
//
- //set daughter
+ //set track -
//
fParamM = paramm;
-
}
+void AliESDV0MI::SetRp(const Double_t *rp){
+ //
+ // set pid +
+ //
+ for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
+}
+
+void AliESDV0MI::SetRm(const Double_t *rm){
+ //
+ // set pid -
+ //
+ for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
+}
+
+
void AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
{
//
//
// updates Kink Info
//
- Float_t distance1,distance2;
+ // Float_t distance1,distance2;
+ Float_t distance2;
//
AliHelix phelix(fParamP);
AliHelix mhelix(fParamM);
Double_t phase[2][2],radius[2];
Int_t points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
Double_t delta1=10000,delta2=10000;
-
+ /*
if (points<=0) return;
if (points>0){
phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
}
distance1 = TMath::Min(delta1,delta2);
+ */
//
//find intersection parabolic
//
fDist2 = TMath::Sqrt(distance2);
//
//
- Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
- Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
- Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
- Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+ Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
+ Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
+ Double_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+ Double_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
vnorm2 = TMath::Sqrt(vnorm2);
- Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
- Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+ Double_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+ Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
pnorm2 = TMath::Sqrt(pnorm2);
fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
// friend class AliITStrackerMI;
AliESDV0MI(); //constructor
//
+ const AliExternalTrackParam *GetParamP() const {return &fParamP;}
+ const AliExternalTrackParam *GetParamM() const {return &fParamM;}
void SetP(const AliExternalTrackParam & paramp);
void SetM(const AliExternalTrackParam & paramd);
+ void SetRp(const Double_t *rp);
+ void SetRm(const Double_t *rm);
void UpdatePID(Double_t pidp[5], Double_t pidm[5]);
+ void SetStatus(Int_t status){fStatus=status;}
+ Int_t GetStatus() const {return fStatus;}
Float_t GetEffMass(UInt_t p1, UInt_t p2);
Float_t GetProb(UInt_t p1, UInt_t p2);
void Update(Float_t vertex[3]); //update
Float_t GetNBefore() const {return fNBefore;}
void SetNBefore(Float_t nb) {fNBefore=nb;}
void SetLab(Int_t i, Int_t lab) {fLab[i]=lab;}
-
+ void SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1);
+ const Float_t * GetCausalityP() const {return fCausality;}
private:
AliExternalTrackParam fParamP;
AliExternalTrackParam fParamM;
Double_t fXr[3]; //rec. position according helix
Double_t fAngle[3]; //three angles
Double_t fRr; //rec position of the vertex
- Int_t fStatus; //status
+ Int_t fStatus; //status - 1 - TPC V0 2- ITS V0 4- accepted - 0 -rejected
Int_t fRow0; // critical layer
Int_t fOrder[3]; //order of the vertex
// quality information
Double_t fDistNorm; //normalized DCA
Double_t fDistSigma; //sigma of distance
+ Float_t fCausality[4]; // causality information - see comments in SetCausality
Float_t fChi2Before; //chi2 of the tracks before V0
Float_t fNBefore; // number of possible points before V0
Float_t fChi2After; // chi2 of the tracks after V0
Float_t fPointAngleTh; //point angle theta
Float_t fPointAngle; //point angle full
-
-
- ClassDef(AliESDV0MI,1) // ESD V0 vertex
+ ClassDef(AliESDV0MI,2) // ESD V0 vertex
};
return fTPCncls;
}
+Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
+ //
+ // GetDensity of the clusters on given region between row0 and row1
+ // Dead zone effect takin into acoount
+ //
+ Int_t good = 0;
+ Int_t found = 0;
+ //
+ for (Int_t i=row0;i<=row1;i++){
+ Int_t index = fTPCindex[i];
+ if (index!=-1) good++; // track outside of dead zone
+ if (index>0) found++;
+ }
+ Float_t density=0.5;
+ if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
+ return density;
+}
//_______________________________________________________________________
void AliESDtrack::SetTPCpid(const Double_t *p) {
// Sets values for the probability of each particle type (in TPC)
- for (Int_t i=0; i<AliPID::kSPECIES; i++) fTPCr[i]=p[i];
+ // normalize probabiluty to 1
+ //
+ //
+// Double_t sump=0;
+// for (Int_t i=0; i<AliPID::kSPECIES; i++) {
+// sump+=p[i];
+// }
+// for (Int_t i=0; i<AliPID::kSPECIES; i++) {
+// if (sump>0){
+// fTPCr[i]=p[i]/sump;
+// }
+// else{
+// fTPCr[i]=p[i];
+// }
+// }
+ for (Int_t i=0; i<AliPID::kSPECIES; i++) fTPCr[i] = p[i];
SetStatus(AliESDtrack::kTPCpid);
}
Float_t GetTPCsignal() const {return fTPCsignal;}
Float_t GetTPCchi2() const {return fTPCchi2;}
Int_t GetTPCclusters(Int_t *idx) const;
+ Float_t GetTPCdensity(Int_t row0, Int_t row1) const;
Int_t GetTPCLabel() const {return fTPCLabel;}
Int_t GetKinkIndex(Int_t i) const { return fKinkIndexes[i];}
Int_t GetV0Index(Int_t i) const { return fV0Indexes[i];}
};
protected:
- AliESDtrack & operator=(const AliESDtrack & );
+ //AliESDtrack & operator=(const AliESDtrack & );
ULong_t fFlags; // Reconstruction status flags
Int_t fLabel; // Track label
// TPC related track information
Float_t fTPCchi2; // chi2 in the TPC
Int_t fTPCncls; // number of clusters assigned in the TPC
- Int_t fTPCindex[180]; //! indices of the assigned TPC clusters
+ Int_t fTPCindex[180]; // indices of the assigned TPC clusters
TBits fTPCClusterMap; // Map of clusters, one bit per padrow; 1 if has a cluster on given padrow
Float_t fTPCsignal; // detector's PID signal
Float_t fTPCr[AliPID::kSPECIES]; // "detector response probabilities" (for the PID)
Float_t fTOFsignal; // detector's PID signal
Float_t fTOFr[AliPID::kSPECIES]; // "detector response probabilities" (for the PID)
Int_t fTOFLabel[3]; // TOF label
- Float_t fTOFInfo[10]; //! TOF informations
+ Float_t fTOFInfo[10]; // TOF informations
// PHOS related track information
Float_t fPHOSpos[3]; // position localised by PHOS in global coordinate system
Float_t fRICHdx; // x of the track impact minus x of the MIP
Float_t fRICHdy; // y of the track impact minus y of the MIP
- ClassDef(AliESDtrack,13) //ESDtrack
+ ClassDef(AliESDtrack,14) //ESDtrack
};
#endif
Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
Int_t GetNindex() const {return fNidx;}
Int_t GetPindex() const {return fPidx;}
+ void SetESDindexes(Int_t ip, Int_t im){fNidx=ip;fPidx=im;}
void SetDcaDaughters(Double_t rDcaDaughters=0.);
Double_t GetDcaDaughters() {return fDcaDaughters;}
-
protected:
Int_t fPdgCode; // reconstructed V0's type (PDG code)
Double_t fEffMass; // reconstructed V0's effective mass
fTRDReferences = new TClonesArray("AliTrackReference",10);
fTOFReferences = new TClonesArray("AliTrackReference",10);
fTRdecay.SetTrack(-1);
+ fCharge = 0;
}
AliMCInfo::~AliMCInfo()
fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
fPointAngle = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
- Float_t mass1 = fMCd.fMass;
- Float_t mass2 = fMCm.fMass;
+ Double_t mass1 = fMCd.fMass;
+ Double_t mass2 = fMCm.fMass;
- Float_t e1 = TMath::Sqrt(mass1*mass1+
- fMCPdr[0]*fMCPdr[0]+
- fMCPdr[1]*fMCPdr[1]+
- fMCPdr[2]*fMCPdr[2]);
- Float_t e2 = TMath::Sqrt(mass2*mass2+
+ Double_t e1 = TMath::Sqrt(mass1*mass1+
+ fMCPd[0]*fMCPd[0]+
+ fMCPd[1]*fMCPd[1]+
+ fMCPd[2]*fMCPd[2]);
+ Double_t e2 = TMath::Sqrt(mass2*mass2+
fMCPm[0]*fMCPm[0]+
fMCPm[1]*fMCPm[1]+
fMCPm[2]*fMCPm[2]);
fInvMass =
- (fMCPm[0]+fMCPdr[0])*(fMCPm[0]+fMCPdr[0])+
- (fMCPm[1]+fMCPdr[1])*(fMCPm[1]+fMCPdr[1])+
- (fMCPm[2]+fMCPdr[2])*(fMCPm[2]+fMCPdr[2]);
+ (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
+ (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
+ (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
// fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
fInvMass = (e1+e2)*(e1+e2)-fInvMass;
-
-
/////////////////////////////////////////////////////////////////////////////////
void AliGenKinkInfo::Update()
{
if (dinfo){
v0info->fMCm = (*info);
v0info->fMCd = (*dinfo);
+ v0info->fMotherP = (*motherparticle);
v0info->Update(fVPrim);
br->SetAddress(&v0info);
fTreeV0->Fill();
public:
AliMCInfo fMCd; //info about daughter particle - second particle for V0
AliMCInfo fMCm; //info about mother particle - first particle for V0
+ TParticle fMotherP; //particle info about mother particle
void Update(Float_t vertex[3]); // put some derived info to special field
Double_t fMCDist1; //info about closest distance according closest MC - linear DCA
Double_t fMCDist2; //info about closest distance parabolic DCA
}
-void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion)
+void AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr)
{
// return momentum at given phase
Double_t x[3],g[3],gg[3];
p[0] = fHelix[8]*g[0]/(mt*conversion);
p[1] = fHelix[8]*g[1]/(mt*conversion);
p[2] = fHelix[8]*g[2]/(mt*conversion);
+ if (xr){
+ xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2];
+ }
}
void AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3])
return 1;
}
+Double_t AliHelix::GetPointAngle(AliHelix &h, Double_t phase[2], const Float_t * vertex)
+{
+ //
+ // get point angle bettwen two helixes
+ //
+ Double_t r0[3],p0[4];
+ Double_t r1[3],p1[4];
+ GetMomentum(phase[0],p0,1,r0);
+ h.GetMomentum(phase[1],p1,1,r1);
+ //
+ Double_t r[3] = {(r0[0]+r1[0])*0.5-vertex[0],(r0[1]+r1[1])*0.5-vertex[1],(r0[2]+r1[2])*0.5-vertex[2]};
+ //intersection point - relative to the prim vertex
+ Double_t p[3] = { p0[0]+p1[0], p0[1]+p1[1],p0[2]+p1[2]};
+ // derivation vector
+ Double_t normr = TMath::Sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);
+ Double_t normp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+ Double_t pointAngle = (r[0]*p[0]+r[1]*p[1]+r[2]*p[2])/(normr*normp);
+ return pointAngle;
+}
+
Double_t AliHelix::GetPhase(Double_t x, Double_t y )
{
void Evaluate(Double_t t, Double_t r[3], //radius vector
Double_t g[3], //first defivatives
Double_t gg[3]); //second derivatives
- void GetMomentum(Double_t phase, Double_t p[4], Double_t conversion=0.); // return momentum
+ void GetMomentum(Double_t phase, Double_t p[4], Double_t conversion=0., Double_t *xr=0); // return momentum
void GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3]);
inline Double_t GetHelixR(Double_t phase=0);
inline Double_t GetHelixZ(Double_t phase=0);
Int_t GetPhase(Double_t r0, Double_t t[2]); //return phase for the nearest point
Int_t GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut=3.);
Int_t GetClosestPhases(AliHelix &h, Double_t phase[2][2]);
+ Double_t GetPointAngle(AliHelix &h, Double_t phase[2],const Float_t *vertex);
Int_t LinearDCA(AliHelix &h, Double_t &t1, Double_t &t2,
Double_t &R, Double_t &dist);
//
//
// Constant L3 field , if this option was selected
//
- b[2] = - fSolenoid;
+ b[2] = (- fSolenoid)*fFactor;
return;
}
} else if (fFieldMap[1]->Inside(xm[0], xm[1], xm[2])) {
TTreeDataElement:: TTreeDataElement(TDataType* type) :
fName(),
- fType(' '),
+ fType(0),
fDType(type),
fClass(0),
fPointer(0)
TTreeDataElement:: TTreeDataElement(TClass* cl) :
fName(),
- fType(' '),
+ fType(0),
fDType(0),
fClass(cl),
fPointer(0)
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+
+
+
+//-----------------------------------------------------------------
+// Implementation of the TPC seed class
+// This class is used by the AliTPCtrackerMI class
+// Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
+//-----------------------------------------------------------------
+#include "TClonesArray.h"
+#include "AliTPCseed.h"
+
+ClassImp(AliTPCseed)
+
+
+
+AliTPCseed::AliTPCseed():AliTPCtrack(){
+ //
+ fRow=0;
+ fRemoval =0;
+ for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
+ for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+ for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
+ for (Int_t i=0;i<5;i++) fTPCr[i]=0.2;
+ fPoints = 0;
+ fEPoints = 0;
+ fNFoundable =0;
+ fNShared =0;
+ fRemoval = 0;
+ fSort =0;
+ fFirstPoint =0;
+ fNoCluster =0;
+ fBSigned = kFALSE;
+ fSeed1 =-1;
+ fSeed2 =-1;
+ fCurrentCluster =0;
+ fCurrentSigmaY2=0;
+ fCurrentSigmaZ2=0;
+ fCircular = 0; // not curling track
+}
+AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
+ //---------------------
+ // dummy copy constructor
+ //-------------------------
+ for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i];
+ for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
+
+ fPoints = 0;
+ fEPoints = 0;
+ fCircular =0;
+}
+AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
+ //
+ //copy constructor
+ fPoints = 0;
+ fEPoints = 0;
+ fNShared =0;
+ // fTrackPoints =0;
+ fRemoval =0;
+ fSort =0;
+ for (Int_t i=0;i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i);
+
+ for (Int_t i=0;i<160;i++) {
+ fClusterPointer[i] = 0;
+ Int_t index = t.GetClusterIndex(i);
+ if (index>=-1){
+ SetClusterIndex2(i,index);
+ }
+ else{
+ SetClusterIndex2(i,-3);
+ }
+ }
+ fFirstPoint =0;
+ fNoCluster =0;
+ fBSigned = kFALSE;
+ fSeed1 =-1;
+ fSeed2 =-1;
+ fCurrentCluster =0;
+ fCurrentSigmaY2=0;
+ fCurrentSigmaZ2=0;
+ fCircular =0;
+}
+
+AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
+ Double_t xr, Double_t alpha):
+ AliTPCtrack(index, xx, cc, xr, alpha) {
+ //
+ //
+ //constructor
+ fRow =0;
+ for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
+ for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+ for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
+ for (Int_t i=0;i<5;i++) fTPCr[i]=0.2;
+
+ fPoints = 0;
+ fEPoints = 0;
+ fNFoundable =0;
+ fNShared = 0;
+ // fTrackPoints =0;
+ fRemoval =0;
+ fSort =0;
+ fFirstPoint =0;
+ // fHelixIn = new TClonesArray("AliHelix",0);
+ //fHelixOut = new TClonesArray("AliHelix",0);
+ fNoCluster =0;
+ fBSigned = kFALSE;
+ fSeed1 =-1;
+ fSeed2 =-1;
+ fCurrentCluster =0;
+ fCurrentSigmaY2=0;
+ fCurrentSigmaZ2=0;
+}
+
+AliTPCseed::~AliTPCseed(){
+ //
+ // destructor
+ if (fPoints) delete fPoints;
+ fPoints =0;
+ if (fEPoints) delete fEPoints;
+ fEPoints = 0;
+ fNoCluster =0;
+}
+
+AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
+{
+ //
+ //
+ return &fTrackPoints[i];
+}
+
+void AliTPCseed::RebuildSeed()
+{
+ //
+ // rebuild seed to be ready for storing
+ AliTPCclusterMI cldummy;
+ cldummy.SetQ(0);
+ AliTPCTrackPoint pdummy;
+ pdummy.GetTPoint().fIsShared = 10;
+ for (Int_t i=0;i<160;i++){
+ AliTPCclusterMI * cl0 = fClusterPointer[i];
+ AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
+ if (cl0){
+ trpoint->GetTPoint() = *(GetTrackPoint(i));
+ trpoint->GetCPoint() = *cl0;
+ trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
+ }
+ else{
+ *trpoint = pdummy;
+ trpoint->GetCPoint()= cldummy;
+ }
+
+ }
+
+}
+
+
+Double_t AliTPCseed::GetDensityFirst(Int_t n)
+{
+ //
+ //
+ // return cluster for n rows bellow first point
+ Int_t nfoundable = 1;
+ Int_t nfound = 1;
+ for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
+ Int_t index = GetClusterIndex2(i);
+ if (index!=-1) nfoundable++;
+ if (index>0) nfound++;
+ }
+ if (nfoundable<n) return 0;
+ return Double_t(nfound)/Double_t(nfoundable);
+
+}
+
+
+void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
+{
+ // get cluster stat. on given region
+ //
+ found = 0;
+ foundable = 0;
+ shared =0;
+ for (Int_t i=first;i<last; i++){
+ Int_t index = GetClusterIndex2(i);
+ if (index!=-1) foundable++;
+ if (fClusterPointer[i]) {
+ found++;
+ }
+ else
+ continue;
+
+ if (fClusterPointer[i]->IsUsed(10)) {
+ shared++;
+ continue;
+ }
+ if (!plus2) continue; //take also neighborhoud
+ //
+ if ( (i>0) && fClusterPointer[i-1]){
+ if (fClusterPointer[i-1]->IsUsed(10)) {
+ shared++;
+ continue;
+ }
+ }
+ if ( fClusterPointer[i+1]){
+ if (fClusterPointer[i+1]->IsUsed(10)) {
+ shared++;
+ continue;
+ }
+ }
+
+ }
+ //if (shared>found){
+ //Error("AliTPCseed::GetClusterStatistic","problem\n");
+ //}
+}
+
+
+
+
+
+void AliTPCseed::Reset(Bool_t all)
+{
+ //
+ //
+ SetNumberOfClusters(0);
+ fNFoundable = 0;
+ SetChi2(0);
+ ResetCovariance();
+ /*
+ if (fTrackPoints){
+ for (Int_t i=0;i<8;i++){
+ delete [] fTrackPoints[i];
+ }
+ delete fTrackPoints;
+ fTrackPoints =0;
+ }
+ */
+
+ if (all){
+ for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
+ for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+ }
+
+}
+
+
+void AliTPCseed::Modify(Double_t factor)
+{
+
+ //------------------------------------------------------------------
+ //This function makes a track forget its history :)
+ //------------------------------------------------------------------
+ if (factor<=0) {
+ ResetCovariance();
+ return;
+ }
+ fC00*=factor;
+ fC10*=0; fC11*=factor;
+ fC20*=0; fC21*=0; fC22*=factor;
+ fC30*=0; fC31*=0; fC32*=0; fC33*=factor;
+ fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor;
+ SetNumberOfClusters(0);
+ fNFoundable =0;
+ SetChi2(0);
+ fRemoval = 0;
+ fCurrentSigmaY2 = 0.000005;
+ fCurrentSigmaZ2 = 0.000005;
+ fNoCluster = 0;
+ //fFirstPoint = 160;
+ //fLastPoint = 0;
+}
+
+
+
+
+Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
+{
+ //-----------------------------------------------------------------
+ // This function find proloncation of a track to a reference plane x=xk.
+ // doesn't change internal state of the track
+ //-----------------------------------------------------------------
+
+ Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
+
+ if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
+ return 0;
+ }
+
+ // Double_t y1=fP0, z1=fP1;
+ Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
+ Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
+
+ y = fP0;
+ z = fP1;
+ //y += dx*(c1+c2)/(r1+r2);
+ //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
+
+ Double_t dy = dx*(c1+c2)/(r1+r2);
+ Double_t dz = 0;
+ //
+ Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
+ /*
+ if (TMath::Abs(delta)>0.0001){
+ dz = fP3*TMath::ASin(delta)/fP4;
+ }else{
+ dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
+ }
+ */
+ // dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4;
+ dz = fP3*TMath::ASin(delta)/fP4;
+ //
+ y+=dy;
+ z+=dz;
+
+
+ return 1;
+}
+
+
+//_____________________________________________________________________________
+Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
+{
+ //-----------------------------------------------------------------
+ // This function calculates a predicted chi2 increment.
+ //-----------------------------------------------------------------
+ //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+ Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
+ r00+=fC00; r01+=fC10; r11+=fC11;
+
+ Double_t det=r00*r11 - r01*r01;
+ if (TMath::Abs(det) < 1.e-10) {
+ //Int_t n=GetNumberOfClusters();
+ //if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+ return 1e10;
+ }
+ Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
+
+ Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+
+ return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
+}
+
+
+//_________________________________________________________________________________________
+
+
+Int_t AliTPCseed::Compare(const TObject *o) const {
+ //-----------------------------------------------------------------
+ // This function compares tracks according to the sector - for given sector according z
+ //-----------------------------------------------------------------
+ AliTPCseed *t=(AliTPCseed*)o;
+
+ if (fSort == 0){
+ if (t->fRelativeSector>fRelativeSector) return -1;
+ if (t->fRelativeSector<fRelativeSector) return 1;
+ Double_t z2 = t->GetZ();
+ Double_t z1 = GetZ();
+ if (z2>z1) return 1;
+ if (z2<z1) return -1;
+ return 0;
+ }
+ else {
+ Float_t f2 =1;
+ f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
+ if (t->fBConstrain) f2=1.2;
+
+ Float_t f1 =1;
+ f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
+
+ if (fBConstrain) f1=1.2;
+
+ if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
+ else return +1;
+ }
+}
+
+
+
+
+//_____________________________________________________________________________
+Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
+ //-----------------------------------------------------------------
+ // This function associates a cluster with this track.
+ //-----------------------------------------------------------------
+ Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
+
+ r00+=fC00; r01+=fC10; r11+=fC11;
+ Double_t det=r00*r11 - r01*r01;
+ Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+
+ Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
+ Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
+ Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
+ Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
+ Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
+
+ Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+ Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
+ if (TMath::Abs(cur*fX-eta) >= 0.9) {
+ return 0;
+ }
+
+ fP0 += k00*dy + k01*dz;
+ fP1 += k10*dy + k11*dz;
+ fP2 = eta;
+ fP3 += k30*dy + k31*dz;
+ fP4 = cur;
+
+ Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
+ Double_t c12=fC21, c13=fC31, c14=fC41;
+
+ fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
+ fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
+ fC40-=k00*c04+k01*c14;
+
+ fC11-=k10*c01+k11*fC11;
+ fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
+ fC41-=k10*c04+k11*c14;
+
+ fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
+ fC42-=k20*c04+k21*c14;
+
+ fC33-=k30*c03+k31*c13;
+ fC43-=k40*c03+k41*c13;
+
+ fC44-=k40*c04+k41*c14;
+
+ Int_t n=GetNumberOfClusters();
+ // fIndex[n]=index;
+ SetNumberOfClusters(n+1);
+ SetChi2(GetChi2()+chisq);
+
+ return 1;
+}
+
+
+
+//_____________________________________________________________________________
+void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
+ //-----------------------------------------------------------------
+ // This funtion calculates dE/dX within the "low" and "up" cuts.
+ //-----------------------------------------------------------------
+
+ Float_t amp[200];
+ Float_t angular[200];
+ Float_t weight[200];
+ Int_t index[200];
+ //Int_t nc = 0;
+ // TClonesArray & arr = *fPoints;
+ Float_t meanlog = 100.;
+
+ Float_t mean[4] = {0,0,0,0};
+ Float_t sigma[4] = {1000,1000,1000,1000};
+ Int_t nc[4] = {0,0,0,0};
+ Float_t norm[4] = {1000,1000,1000,1000};
+ //
+ //
+ fNShared =0;
+
+ for (Int_t of =0; of<4; of++){
+ for (Int_t i=of+i1;i<i2;i+=4)
+ {
+ Int_t index = fIndex[i];
+ if (index<0||index&0x8000) continue;
+
+ //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
+ AliTPCTrackerPoint * point = GetTrackPoint(i);
+ //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
+ //AliTPCTrackerPoint * pointp = 0;
+ //if (i<159) pointp = GetTrackPoint(i+1);
+
+ if (point==0) continue;
+ AliTPCclusterMI * cl = fClusterPointer[i];
+ if (cl==0) continue;
+ if (onlyused && (!cl->IsUsed(10))) continue;
+ if (cl->IsUsed(11)) {
+ fNShared++;
+ continue;
+ }
+ Int_t type = cl->GetType();
+ //if (point->fIsShared){
+ // fNShared++;
+ // continue;
+ //}
+ //if (pointm)
+ // if (pointm->fIsShared) continue;
+ //if (pointp)
+ // if (pointp->fIsShared) continue;
+
+ if (type<0) continue;
+ //if (type>10) continue;
+ //if (point->GetErrY()==0) continue;
+ //if (point->GetErrZ()==0) continue;
+
+ //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
+ //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
+ //if ((ddy*ddy+ddz*ddz)>10) continue;
+
+
+ // if (point->GetCPoint().GetMax()<5) continue;
+ if (cl->GetMax()<5) continue;
+ Float_t angley = point->GetAngleY();
+ Float_t anglez = point->GetAngleZ();
+
+ Float_t rsigmay2 = point->GetSigmaY();
+ Float_t rsigmaz2 = point->GetSigmaZ();
+ /*
+ Float_t ns = 1.;
+ if (pointm){
+ rsigmay += pointm->GetTPoint().GetSigmaY();
+ rsigmaz += pointm->GetTPoint().GetSigmaZ();
+ ns+=1.;
+ }
+ if (pointp){
+ rsigmay += pointp->GetTPoint().GetSigmaY();
+ rsigmaz += pointp->GetTPoint().GetSigmaZ();
+ ns+=1.;
+ }
+ rsigmay/=ns;
+ rsigmaz/=ns;
+ */
+
+ Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
+
+ Float_t ampc = 0; // normalization to the number of electrons
+ if (i>64){
+ // ampc = 1.*point->GetCPoint().GetMax();
+ ampc = 1.*cl->GetMax();
+ //ampc = 1.*point->GetCPoint().GetQ();
+ // AliTPCClusterPoint & p = point->GetCPoint();
+ // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
+ // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
+ //Float_t dz =
+ // TMath::Abs( Int_t(iz) - iz + 0.5);
+ //ampc *= 1.15*(1-0.3*dy);
+ //ampc *= 1.15*(1-0.3*dz);
+ // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
+ //ampc *=zfactor;
+ }
+ else{
+ //ampc = 1.0*point->GetCPoint().GetMax();
+ ampc = 1.0*cl->GetMax();
+ //ampc = 1.0*point->GetCPoint().GetQ();
+ //AliTPCClusterPoint & p = point->GetCPoint();
+ // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
+ //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
+ //Float_t dz =
+ // TMath::Abs( Int_t(iz) - iz + 0.5);
+
+ //ampc *= 1.15*(1-0.3*dy);
+ //ampc *= 1.15*(1-0.3*dz);
+ // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
+ //ampc *=zfactor;
+
+ }
+ ampc *= 2.0; // put mean value to channel 50
+ //ampc *= 0.58; // put mean value to channel 50
+ Float_t w = 1.;
+ // if (type>0) w = 1./(type/2.-0.5);
+ // Float_t z = TMath::Abs(cl->GetZ());
+ if (i<64) {
+ ampc /= 0.6;
+ //ampc /= (1+0.0008*z);
+ } else
+ if (i>128){
+ ampc /=1.5;
+ //ampc /= (1+0.0008*z);
+ }else{
+ //ampc /= (1+0.0008*z);
+ }
+
+ if (type<0) { //amp at the border - lower weight
+ // w*= 2.;
+
+ continue;
+ }
+ if (rsigma>1.5) ampc/=1.3; // if big backround
+ amp[nc[of]] = ampc;
+ angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
+ weight[nc[of]] = w;
+ nc[of]++;
+ }
+
+ TMath::Sort(nc[of],amp,index,kFALSE);
+ Float_t sumamp=0;
+ Float_t sumamp2=0;
+ Float_t sumw=0;
+ //meanlog = amp[index[Int_t(nc[of]*0.33)]];
+ meanlog = 50;
+ for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
+ Float_t ampl = amp[index[i]]/angular[index[i]];
+ ampl = meanlog*TMath::Log(1.+ampl/meanlog);
+ //
+ sumw += weight[index[i]];
+ sumamp += weight[index[i]]*ampl;
+ sumamp2 += weight[index[i]]*ampl*ampl;
+ norm[of] += angular[index[i]]*weight[index[i]];
+ }
+ if (sumw<1){
+ SetdEdx(0);
+ }
+ else {
+ norm[of] /= sumw;
+ mean[of] = sumamp/sumw;
+ sigma[of] = sumamp2/sumw-mean[of]*mean[of];
+ if (sigma[of]>0.1)
+ sigma[of] = TMath::Sqrt(sigma[of]);
+ else
+ sigma[of] = 1000;
+
+ mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
+ //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
+ //mean *=(1-0.1*(norm-1.));
+ }
+ }
+
+ Float_t dedx =0;
+ fSdEdx =0;
+ fMAngular =0;
+ // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
+ // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
+
+
+ // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
+ // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
+
+ Int_t norm2 = 0;
+ Int_t norm3 = 0;
+ for (Int_t i =0;i<4;i++){
+ if (nc[i]>2&&nc[i]<1000){
+ dedx += mean[i] *nc[i];
+ fSdEdx += sigma[i]*(nc[i]-2);
+ fMAngular += norm[i] *nc[i];
+ norm2 += nc[i];
+ norm3 += nc[i]-2;
+ }
+ fDEDX[i] = mean[i];
+ fSDEDX[i] = sigma[i];
+ fNCDEDX[i]= nc[i];
+ }
+
+ if (norm3>0){
+ dedx /=norm2;
+ fSdEdx /=norm3;
+ fMAngular/=norm2;
+ }
+ else{
+ SetdEdx(0);
+ return;
+ }
+ // Float_t dedx1 =dedx;
+ /*
+ dedx =0;
+ for (Int_t i =0;i<4;i++){
+ if (nc[i]>2&&nc[i]<1000){
+ mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
+ dedx += mean[i] *nc[i];
+ }
+ fDEDX[i] = mean[i];
+ }
+ dedx /= norm2;
+ */
+
+
+ SetdEdx(dedx);
+
+ //mi deDX
+
+
+
+ //Very rough PID
+ Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
+
+ if (p<0.6) {
+ if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+ if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
+ SetMass(0.93827); return;
+ }
+
+ if (p<1.2) {
+ if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+ SetMass(0.93827); return;
+ }
+
+ SetMass(0.13957); return;
+
+}
+Double_t AliTPCseed::Bethe(Double_t bg){
+ //
+ // This is the Bethe-Bloch function normalised to 1 at the minimum
+ //
+ Double_t bg2=bg*bg;
+ Double_t bethe;
+ if (bg<3.5e1)
+ bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
+ else // Density effect ( approximately :)
+ bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2));
+ return bethe/11.091;
+}
+
+void AliTPCseed::CookPID()
+{
+ //
+ // cook PID information according dEdx
+ //
+ Double_t fRange = 10.;
+ Double_t fRes = 0.1;
+ Double_t fMIP = 47.;
+ //
+ Int_t ns=AliPID::kSPECIES;
+ Double_t sumr =0;
+ for (Int_t j=0; j<ns; j++) {
+ Double_t mass=AliPID::ParticleMass(j);
+ Double_t mom=P();
+ Double_t dedx=fdEdx/fMIP;
+ Double_t bethe=Bethe(mom/mass);
+ Double_t sigma=fRes*bethe;
+ if (sigma>0.001){
+ if (TMath::Abs(dedx-bethe) > fRange*sigma) {
+ fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
+ sumr+=fTPCr[j];
+ continue;
+ }
+ fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
+ sumr+=fTPCr[j];
+ }
+ else{
+ fTPCr[j]=1.;
+ sumr+=fTPCr[j];
+ }
+ }
+ for (Int_t j=0; j<ns; j++) {
+ fTPCr[j]/=sumr; //normalize
+ }
+}
+
+/*
+void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
+ //-----------------------------------------------------------------
+ // This funtion calculates dE/dX within the "low" and "up" cuts.
+ //-----------------------------------------------------------------
+
+ Float_t amp[200];
+ Float_t angular[200];
+ Float_t weight[200];
+ Int_t index[200];
+ Bool_t inlimit[200];
+ for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
+ for (Int_t i=0;i<200;i++) amp[i]=10000;
+ for (Int_t i=0;i<200;i++) angular[i]= 1;;
+
+
+ //
+ Float_t meanlog = 100.;
+ Int_t indexde[4]={0,64,128,160};
+
+ Float_t amean =0;
+ Float_t asigma =0;
+ Float_t anc =0;
+ Float_t anorm =0;
+
+ Float_t mean[4] = {0,0,0,0};
+ Float_t sigma[4] = {1000,1000,1000,1000};
+ Int_t nc[4] = {0,0,0,0};
+ Float_t norm[4] = {1000,1000,1000,1000};
+ //
+ //
+ fNShared =0;
+
+ // for (Int_t of =0; of<3; of++){
+ // for (Int_t i=indexde[of];i<indexde[of+1];i++)
+ for (Int_t i =0; i<160;i++)
+ {
+ AliTPCTrackPoint * point = GetTrackPoint(i);
+ if (point==0) continue;
+ if (point->fIsShared){
+ fNShared++;
+ continue;
+ }
+ Int_t type = point->GetCPoint().GetType();
+ if (type<0) continue;
+ if (point->GetCPoint().GetMax()<5) continue;
+ Float_t angley = point->GetTPoint().GetAngleY();
+ Float_t anglez = point->GetTPoint().GetAngleZ();
+ Float_t rsigmay = point->GetCPoint().GetSigmaY();
+ Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
+ Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
+
+ Float_t ampc = 0; // normalization to the number of electrons
+ if (i>64){
+ ampc = point->GetCPoint().GetMax();
+ }
+ else{
+ ampc = point->GetCPoint().GetMax();
+ }
+ ampc *= 2.0; // put mean value to channel 50
+ // ampc *= 0.565; // put mean value to channel 50
+
+ Float_t w = 1.;
+ Float_t z = TMath::Abs(point->GetCPoint().GetZ());
+ if (i<64) {
+ ampc /= 0.63;
+ } else
+ if (i>128){
+ ampc /=1.51;
+ }
+ if (type<0) { //amp at the border - lower weight
+ continue;
+ }
+ if (rsigma>1.5) ampc/=1.3; // if big backround
+ angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
+ amp[i] = ampc/angular[i];
+ weight[i] = w;
+ anc++;
+ }
+
+ TMath::Sort(159,amp,index,kFALSE);
+ for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
+ inlimit[index[i]] = kTRUE; // take all clusters
+ }
+
+ // meanlog = amp[index[Int_t(anc*0.3)]];
+ meanlog =10000.;
+ for (Int_t of =0; of<3; of++){
+ Float_t sumamp=0;
+ Float_t sumamp2=0;
+ Float_t sumw=0;
+ for (Int_t i=indexde[of];i<indexde[of+1];i++)
+ {
+ if (inlimit[i]==kFALSE) continue;
+ Float_t ampl = amp[i];
+ ///angular[i];
+ ampl = meanlog*TMath::Log(1.+ampl/meanlog);
+ //
+ sumw += weight[i];
+ sumamp += weight[i]*ampl;
+ sumamp2 += weight[i]*ampl*ampl;
+ norm[of] += angular[i]*weight[i];
+ nc[of]++;
+ }
+ if (sumw<1){
+ SetdEdx(0);
+ }
+ else {
+ norm[of] /= sumw;
+ mean[of] = sumamp/sumw;
+ sigma[of] = sumamp2/sumw-mean[of]*mean[of];
+ if (sigma[of]>0.1)
+ sigma[of] = TMath::Sqrt(sigma[of]);
+ else
+ sigma[of] = 1000;
+ mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
+ }
+ }
+
+ Float_t dedx =0;
+ fSdEdx =0;
+ fMAngular =0;
+ //
+ Int_t norm2 = 0;
+ Int_t norm3 = 0;
+ Float_t www[3] = {12.,14.,17.};
+ //Float_t www[3] = {1.,1.,1.};
+
+ for (Int_t i =0;i<3;i++){
+ if (nc[i]>2&&nc[i]<1000){
+ dedx += mean[i] *nc[i]*www[i]/sigma[i];
+ fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
+ fMAngular += norm[i] *nc[i];
+ norm2 += nc[i]*www[i]/sigma[i];
+ norm3 += (nc[i]-2)*www[i]/sigma[i];
+ }
+ fDEDX[i] = mean[i];
+ fSDEDX[i] = sigma[i];
+ fNCDEDX[i]= nc[i];
+ }
+
+ if (norm3>0){
+ dedx /=norm2;
+ fSdEdx /=norm3;
+ fMAngular/=norm2;
+ }
+ else{
+ SetdEdx(0);
+ return;
+ }
+ // Float_t dedx1 =dedx;
+
+ dedx =0;
+ Float_t norm4 = 0;
+ for (Int_t i =0;i<3;i++){
+ if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
+ //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
+ dedx += mean[i] *(nc[i])/(sigma[i]);
+ norm4 += (nc[i])/(sigma[i]);
+ }
+ fDEDX[i] = mean[i];
+ }
+ if (norm4>0) dedx /= norm4;
+
+
+
+ SetdEdx(dedx);
+
+ //mi deDX
+
+}
+*/
--- /dev/null
+#ifndef ALITPCSEED_H
+#define ALITPCSEED_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+
+/* $Id$ */
+
+//-------------------------------------------------------
+// TPC seed
+// Class needed for TPC parallel tracking
+//
+// Origin:
+//-------------------------------------------------------
+
+#include <TError.h>
+
+#include "AliTPCtrack.h"
+#include "AliComplexCluster.h"
+#include "AliPID.h"
+
+class TFile;
+class AliTPCParam;
+class AliTPCseed;
+class AliTPCclusterMI;
+class AliTPCTrackerPoint;
+class AliESD;
+class TClonesArray;
+
+class AliTPCseed : public AliTPCtrack {
+ friend class AliTPCtrackerMI;
+ public:
+ AliTPCseed();
+ virtual ~AliTPCseed();
+ AliTPCseed(const AliTPCtrack &t);
+ AliTPCseed(const AliTPCseed &s);
+ //AliTPCseed(const AliTPCseed &t, Double_t a);
+ AliTPCseed(UInt_t index, const Double_t xx[5],
+ const Double_t cc[15], Double_t xr, Double_t alpha);
+ Int_t Compare(const TObject *o) const;
+ void Reset(Bool_t all = kTRUE);
+ Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const;
+ virtual Double_t GetPredictedChi2(const AliTPCclusterMI *cluster) const;
+ virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i);
+ AliTPCTrackerPoint * GetTrackPoint(Int_t i);
+ void RebuildSeed(); // rebuild seed to be ready for storing
+ Double_t GetDensityFirst(Int_t n);
+ Double_t GetSigma2C() const {return fC44;};
+ void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2);
+
+ void Modify(Double_t factor);
+ void SetClusterIndex2(Int_t row, Int_t index) {
+ fIndex[row] = index;
+ }
+ Int_t GetClusterIndex2(Int_t row) const {
+ return fIndex[row];
+ }
+ Int_t GetClusterSector(Int_t row) const {
+ Int_t pica = -1;
+ if (fIndex[row]>=0) pica = ((fIndex[row]&0xff000000)>>24);
+ return pica;
+ }
+
+ void SetErrorY2(Float_t sy2){fErrorY2=sy2;}
+ void SetErrorZ2(Float_t sz2){fErrorZ2=sz2;}
+ void CookdEdx(Double_t low=0.05, Double_t up=0.70, Int_t i1=0, Int_t i2=159, Bool_t onlyused = kFALSE);
+ void CookPID();
+ Double_t Bethe(Double_t bg); // return bethe-bloch
+ // void CookdEdx2(Double_t low=0.05, Double_t up=0.70);
+ Bool_t IsActive() const { return !(fRemoval);}
+ void Desactivate(Int_t reason){ fRemoval = reason;}
+ //
+ //
+ private:
+ // AliTPCseed & operator = (const AliTPCseed &)
+ // {::Fatal("= operator","Not Implemented\n");return *this;}
+ AliESDtrack * fEsd; //!
+ AliTPCclusterMI* fClusterPointer[160]; //! array of cluster pointers -
+ TClonesArray * fPoints; // array with points along the track
+ TClonesArray * fEPoints; // array with exact points - calculated in special macro not used in tracking
+ //---CURRENT VALUES
+ Int_t fRow; //!current row number
+ Int_t fSector; //!current sector number
+ Int_t fRelativeSector; //! index of current relative sector
+ Float_t fCurrentSigmaY2; //!expected current cluster sigma Y
+ Float_t fCurrentSigmaZ2; //!expected current cluster sigma Z
+ Float_t fErrorY2; //!sigma of current cluster
+ Float_t fErrorZ2; //!sigma of current cluster
+ AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation
+ Int_t fCurrentClusterIndex1; //! index of the current cluster
+ Bool_t fInDead; //! indicate if the track is in dead zone
+ Bool_t fIsSeeding; //!indicates if it is proces of seeading
+ Int_t fNoCluster; //!indicates number of rows without clusters
+ Int_t fSort; //!indicate criteria for sorting
+ Bool_t fBSigned; //indicates that clusters of this trackes are signed to be used
+ //
+ //
+ Float_t fDEDX[4]; // dedx according padrows
+ Float_t fSDEDX[4]; // sdedx according padrows
+ Int_t fNCDEDX[4]; // number of clusters for dedx measurment
+ Double_t fTPCr[AliPID::kSPECIES]; // rough PID according TPC
+ //
+ Int_t fSeedType; //seeding type
+ Int_t fSeed1; //first row for seeding
+ Int_t fSeed2; //last row for seeding
+ Int_t fOverlapLabels[12]; //track labels and the length of the overlap
+ Float_t fMAngular; // mean angular factor
+ Char_t fCircular; // indicates curlin track
+ AliTPCTrackerPoint fTrackPoints[160]; //!track points - array track points
+
+ ClassDef(AliTPCseed,1)
+};
+
+
+
+
+
+#endif
+
+
fX = fP0 = fP1 = fP2 = fP3 = fP3 = fP4 = 0.0;
fAlpha = fdEdx = 0.0;
fNumber = 0; // [SR, 01.04.2003]
- for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
+ for (Int_t i=0; i<3;i++) {fKinkIndexes[i]=0; fV0Indexes[i]=0;}
}
//_________________________________________________________________________
fTrackType = 0;
fLab2 = 0;
for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
+ for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
}
//_____________________________________________________________________________
SetLabel(t.GetLabel());
SetMass(t.GetMass());
for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i);
+ for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.GetV0Index(i);
fdEdx = t.GetTPCsignal();
fAlpha = t.GetAlpha();
fTrackType = t.fTrackType;
fLab2 = t.fLab2;
for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i];
+ for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.fV0Indexes[i];
}
//_____________________________________________________________________________
Double_t GetD(Double_t x=0, Double_t y=0) const;
AliExternalTrackParam & GetReference(){ return fReference;}
void UpdateReference(){ new (&fReference) AliExternalTrackParam(*this);}
- Int_t GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];}
+ Int_t GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];}
Int_t* GetKinkIndexes() { return fKinkIndexes;}
+ Int_t GetV0Index(Int_t i) const{ return fV0Indexes[i];}
+ Int_t* GetV0Indexes() { return fV0Indexes;}
protected:
Double_t fX; // X-coordinate of this track (reference plane)
Double_t fAlpha; // Rotation angle the local (TPC sector)
AliExternalTrackParam fReference; // track parameters at the middle of the chamber
Float_t fKinkPoint[12]; //radius, of kink, dfi and dtheta
Int_t fKinkIndexes[3]; // kink indexes - minus = mother + daughter
+ Int_t fV0Indexes[3]; // kink indexes - minus = mother + daughter
ClassDef(AliTPCtrack,2) // Time Projection Chamber reconstructed tracks
};
#include "AliTPCReconstructor.h"
#include "AliTPCclusterMI.h"
#include "AliTPCpolyTrack.h"
-#include "AliTPCreco.h"
+#include "AliTPCreco.h"
+#include "AliTPCseed.h"
#include "AliTPCtrackerMI.h"
#include "TStopwatch.h"
+#include "AliTPCReconstructor.h"
+#include "AliESDkink.h"
+#include "AliPID.h"
#include "TTreeStream.h"
//
-ClassImp(AliTPCseed)
ClassImp(AliTPCtrackerMI)
fNewIO =0;
fDebug =0;
fEvent =0;
+ fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
}
//________________________________________________________________________
AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
fSeeds->Delete();
delete fSeeds;
}
+ if (fDebugStreamer) delete fDebugStreamer;
}
void AliTPCtrackerMI::SetIO()
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+ iotrack.SetV0Indexes(pt->GetV0Indexes());
+ // iotrack.SetTPCpid(pt->fTPCr);
//iotrack.SetTPCindex(i);
fEvent->AddTrack(&iotrack);
continue;
iotrack.SetTPCPoints(pt->GetPoints());
//iotrack.SetTPCindex(i);
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+ iotrack.SetV0Indexes(pt->GetV0Indexes());
+ // iotrack.SetTPCpid(pt->fTPCr);
fEvent->AddTrack(&iotrack);
continue;
}
//iotrack.SetTPCindex(i);
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+ iotrack.SetV0Indexes(pt->GetV0Indexes());
+ //iotrack.SetTPCpid(pt->fTPCr);
fEvent->AddTrack(&iotrack);
continue;
}
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+ iotrack.SetV0Indexes(pt->GetV0Indexes());
+ //iotrack.SetTPCpid(pt->fTPCr);
//iotrack.SetTPCindex(i);
fEvent->AddTrack(&iotrack);
continue;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+ iotrack.SetV0Indexes(pt->GetV0Indexes());
+ //iotrack.SetTPCpid(pt->fTPCr);
//iotrack.SetTPCindex(i);
fEvent->AddTrack(&iotrack);
continue;
iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
iotrack.SetTPCPoints(pt->GetPoints());
iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+ iotrack.SetV0Indexes(pt->GetV0Indexes());
+ // iotrack.SetTPCpid(pt->fTPCr);
//iotrack.SetTPCindex(i);
fEvent->AddTrack(&iotrack);
continue;
-void AliTPCseed::Reset(Bool_t all)
-{
- //
- //
- SetNumberOfClusters(0);
- fNFoundable = 0;
- SetChi2(0);
- ResetCovariance();
- /*
- if (fTrackPoints){
- for (Int_t i=0;i<8;i++){
- delete [] fTrackPoints[i];
- }
- delete fTrackPoints;
- fTrackPoints =0;
- }
- */
-
- if (all){
- for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
- for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
- }
-
-}
-
-
-void AliTPCseed::Modify(Double_t factor)
-{
-
- //------------------------------------------------------------------
- //This function makes a track forget its history :)
- //------------------------------------------------------------------
- if (factor<=0) {
- ResetCovariance();
- return;
- }
- fC00*=factor;
- fC10*=0; fC11*=factor;
- fC20*=0; fC21*=0; fC22*=factor;
- fC30*=0; fC31*=0; fC32*=0; fC33*=factor;
- fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor;
- SetNumberOfClusters(0);
- fNFoundable =0;
- SetChi2(0);
- fRemoval = 0;
- fCurrentSigmaY2 = 0.000005;
- fCurrentSigmaZ2 = 0.000005;
- fNoCluster = 0;
- //fFirstPoint = 160;
- //fLastPoint = 0;
-}
-
-
-
-
-Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
-{
- //-----------------------------------------------------------------
- // This function find proloncation of a track to a reference plane x=xk.
- // doesn't change internal state of the track
- //-----------------------------------------------------------------
-
- Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
-
- if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
- return 0;
- }
-
- // Double_t y1=fP0, z1=fP1;
- Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
- Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
-
- y = fP0;
- z = fP1;
- //y += dx*(c1+c2)/(r1+r2);
- //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
-
- Double_t dy = dx*(c1+c2)/(r1+r2);
- Double_t dz = 0;
- //
- Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
- /*
- if (TMath::Abs(delta)>0.0001){
- dz = fP3*TMath::ASin(delta)/fP4;
- }else{
- dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
- }
- */
- dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4;
- //
- y+=dy;
- z+=dz;
-
-
- return 1;
-}
-
-
-//_____________________________________________________________________________
-Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
-{
- //-----------------------------------------------------------------
- // This function calculates a predicted chi2 increment.
- //-----------------------------------------------------------------
- //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
- Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
- r00+=fC00; r01+=fC10; r11+=fC11;
-
- Double_t det=r00*r11 - r01*r01;
- if (TMath::Abs(det) < 1.e-10) {
- Int_t n=GetNumberOfClusters();
- if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
- return 1e10;
- }
- Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
-
- Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
-
- return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
-}
-
-
-//_________________________________________________________________________________________
-
-
-Int_t AliTPCseed::Compare(const TObject *o) const {
- //-----------------------------------------------------------------
- // This function compares tracks according to the sector - for given sector according z
- //-----------------------------------------------------------------
- AliTPCseed *t=(AliTPCseed*)o;
-
- if (fSort == 0){
- if (t->fRelativeSector>fRelativeSector) return -1;
- if (t->fRelativeSector<fRelativeSector) return 1;
- Double_t z2 = t->GetZ();
- Double_t z1 = GetZ();
- if (z2>z1) return 1;
- if (z2<z1) return -1;
- return 0;
- }
- else {
- Float_t f2 =1;
- f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
- if (t->fBConstrain) f2=1.2;
-
- Float_t f1 =1;
- f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
-
- if (fBConstrain) f1=1.2;
-
- if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
- else return +1;
- }
-}
void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
{
-
-//_____________________________________________________________________________
-Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
- //-----------------------------------------------------------------
- // This function associates a cluster with this track.
- //-----------------------------------------------------------------
- Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
-
- r00+=fC00; r01+=fC10; r11+=fC11;
- Double_t det=r00*r11 - r01*r01;
- Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
-
- Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
- Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
- Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
- Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
- Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
-
- Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
- Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
- if (TMath::Abs(cur*fX-eta) >= 0.9) {
- return 0;
- }
-
- fP0 += k00*dy + k01*dz;
- fP1 += k10*dy + k11*dz;
- fP2 = eta;
- fP3 += k30*dy + k31*dz;
- fP4 = cur;
-
- Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
- Double_t c12=fC21, c13=fC31, c14=fC41;
-
- fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
- fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
- fC40-=k00*c04+k01*c14;
-
- fC11-=k10*c01+k11*fC11;
- fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
- fC41-=k10*c04+k11*c14;
-
- fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
- fC42-=k20*c04+k21*c14;
-
- fC33-=k30*c03+k31*c13;
- fC43-=k40*c03+k41*c13;
-
- fC44-=k40*c04+k41*c14;
-
- Int_t n=GetNumberOfClusters();
- // fIndex[n]=index;
- SetNumberOfClusters(n+1);
- SetChi2(GetChi2()+chisq);
-
- return 1;
-}
-
-
-
//_____________________________________________________________________________
Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
Double_t x2,Double_t y2,
//
TObjArray *kinks= new TObjArray(10000);
+ // TObjArray *v0s= new TObjArray(10000);
Int_t nentries = array->GetEntriesFast();
AliHelix *helixes = new AliHelix[nentries];
Int_t *sign = new Int_t[nentries];
Float_t *fim = new Float_t[nentries];
Float_t *shared = new Float_t[nentries];
Bool_t *circular = new Bool_t[nentries];
+ Float_t *dca = new Float_t[nentries];
+ //const AliESDVertex * primvertex = esd->GetVertex();
//
// nentries = array->GetEntriesFast();
//
usage[i]=0;
AliTPCseed* track = (AliTPCseed*)array->At(i);
if (!track) continue;
+ track->fCircular =0;
shared[i] = kFALSE;
track->UpdatePoints();
if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
}
z0[i]=1000;
circular[i]= kFALSE;
- if (track->GetProlongation(0,y,z)) z0[i] = TMath::Abs(z);
+ if (track->GetProlongation(0,y,z)) z0[i] = z;
+ dca[i] = track->GetD(0,0);
}
//
//
//
// Find circling track
- // TTreeSRedirector cstream("circling.root");
+ TTreeSRedirector &cstream = *fDebugStreamer;
//
for (Int_t i0=0;i0<nentries;i0++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
for (Int_t i1=i0+1;i1<nentries;i1++){
AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
if (!track1) continue;
- if (TMath::Abs(1./track1->fP4)>200) continue;
- if (track1->fN<40) continue;
- if (z0[i0]<20&&z0[i1]<20) continue;
- if (track1->fP4*track0->fP4>0) continue;
- if (track1->fP3*track0->fP3>0) continue;
+ if (track1->fN<40) continue;
+ if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
if (track0->fBConstrain&&track1->fBConstrain) continue;
+ if (TMath::Abs(1./track1->fP4)>200) continue;
+ if (track1->fP4*track0->fP4>0) continue;
+ if (track1->fP3*track0->fP3>0) continue;
if (max(TMath::Abs(1./track0->fP4),TMath::Abs(1./track1->fP4))>190) continue;
if (track0->fBConstrain&&TMath::Abs(track1->fP4)<TMath::Abs(track0->fP4)) continue; //returning - lower momenta
if (track1->fBConstrain&&TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)) continue; //returning - lower momenta
//
- if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
+ Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
+ if (mindcar<5) continue;
+ Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
+ if (mindcaz<5) continue;
+ if (mindcar+mindcaz<20) continue;
+ //
+ //
Float_t xc0 = helixes[i0].GetHelix(6);
Float_t yc0 = helixes[i0].GetHelix(7);
Float_t r0 = helixes[i0].GetHelix(8);
Float_t rmean = (r0+r1)*0.5;
Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
- if (delta>30) continue;
+ //if (delta>30) continue;
if (delta>rmean*0.25) continue;
if (TMath::Abs(r0-r1)/rmean>0.3) continue;
//
*/
if (npoints>0){
Int_t ibest=0;
- helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],5);
+ helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
if (npoints==2){
- helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],5);
+ helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
if (deltah[1]<deltah[0]) ibest=1;
}
deltabest = TMath::Sqrt(deltah[ibest]);
helixes[i0].Evaluate(phase[ibest][0],xyz0);
helixes[i1].Evaluate(phase[ibest][1],xyz1);
helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
- //Double_t radiusbest = TMath::Sqrt(radius[ibest]);
+ Double_t radiusbest = TMath::Sqrt(radius[ibest]);
//
- if (deltabest<25){
- /*
- cstream<<"CC1"<<track0->fLab<<track1->fLab<< //0-1
- track0->fP3<<track1->fP3<< //2-3
- track0->fP4<<track1->fP4<< //4-5
- delta<<rmean<<npoints<< //6-8
- hangles[0]<<hangles[2]<< //9-10
- xyz0[2]<<xyz1[2]<<radiusbest<<deltabest<< //11--14
- track0->fFirstPoint<<track1->fFirstPoint<< //15-16
- track0->fLastPoint<<track1->fLastPoint<< //17-18
- track0<<track1<< //19-20
- phase[ibest][0]<<phase[ibest][1]<<"\n"; //21-22
- */
- }
+ if (deltabest>6) continue;
+ if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
Bool_t sign =kFALSE;
- if (deltabest<3&&hangles[2]>3.06) sign =kTRUE;
- if (TMath::Min(TMath::Abs(track0->GetD()),TMath::Abs(track1->GetD()))>10&&deltabest<6&&hangles[2]>3.) sign= kTRUE;
+ if (hangles[2]>3.06) sign =kTRUE;
+ //
if (sign){
circular[i0] = kTRUE;
- circular[i1] = kTRUE;
+ circular[i1] = kTRUE;
+ if (TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)){
+ track0->fCircular += 1;
+ track1->fCircular += 2;
+ }
+ else{
+ track1->fCircular += 1;
+ track0->fCircular += 2;
+ }
+ }
+ if (sign){
+ //debug stream
+ cstream<<"Curling"<<
+ "lab0="<<track0->fLab<<
+ "lab1="<<track1->fLab<<
+ "Tr0.="<<track0<<
+ "Tr1.="<<track1<<
+ "dca0="<<dca[i0]<<
+ "dca1="<<dca[i1]<<
+ "mindcar="<<mindcar<<
+ "mindcaz="<<mindcaz<<
+ "delta="<<delta<<
+ "rmean="<<rmean<<
+ "npoints="<<npoints<<
+ "hangles0="<<hangles[0]<<
+ "hangles2="<<hangles[2]<<
+ "xyz0="<<xyz0[2]<<
+ "xyzz1="<<xyz1[2]<<
+ "z0="<<z0[i0]<<
+ "z1="<<z0[i1]<<
+ "radius="<<radiusbest<<
+ "deltabest="<<deltabest<<
+ "phase0="<<phase[ibest][0]<<
+ "phase1="<<phase[ibest][1]<<
+ "\n";
}
}
}
}
//
- //
- //
+ // Finf kinks loop
+ //
//
for (Int_t i =0;i<nentries;i++){
if (sign[i]==0) continue;
//
RemoveUsed2(array,0.5,0.4,30);
UnsignClusters();
+ for (Int_t i=0;i<nentries;i++){
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ if (!track0) continue;
+ track0->CookdEdx(0.02,0.6);
+ track0->CookPID();
+ }
//
for (Int_t i=0;i<nentries;i++){
AliTPCseed * track0 = (AliTPCseed*)array->At(i);
delete [] mothers;
//
//
+ delete [] dca;
delete []circular;
delete []shared;
delete []quality;
timer.Print();
}
+void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESD *esd)
+{
+ //
+ // find V0s
+ //
+ //
+ TObjArray *tpcv0s = new TObjArray(100000);
+ Int_t nentries = array->GetEntriesFast();
+ AliHelix *helixes = new AliHelix[nentries];
+ Int_t *sign = new Int_t[nentries];
+ Float_t *alpha = new Float_t[nentries];
+ Float_t *z0 = new Float_t[nentries];
+ Float_t *dca = new Float_t[nentries];
+ Float_t *sdcar = new Float_t[nentries];
+ Float_t *cdcar = new Float_t[nentries];
+ Float_t *pulldcar = new Float_t[nentries];
+ Float_t *pulldcaz = new Float_t[nentries];
+ Float_t *pulldca = new Float_t[nentries];
+ Bool_t *isPrim = new Bool_t[nentries];
+ const AliESDVertex * primvertex = esd->GetVertex();
+ Double_t zvertex = primvertex->GetZv();
+ //
+ // nentries = array->GetEntriesFast();
+ //
+ for (Int_t i=0;i<nentries;i++){
+ sign[i]=0;
+ isPrim[i]=0;
+ AliTPCseed* track = (AliTPCseed*)array->At(i);
+ if (!track) continue;
+ track->GetV0Indexes()[0] = 0; //rest v0 indexes
+ track->GetV0Indexes()[1] = 0; //rest v0 indexes
+ track->GetV0Indexes()[2] = 0; //rest v0 indexes
+ //
+ alpha[i] = track->GetAlpha();
+ new (&helixes[i]) AliHelix(*track);
+ Double_t xyz[3];
+ helixes[i].Evaluate(0,xyz);
+ sign[i] = (track->GetC()>0) ? -1:1;
+ Double_t x,y,z;
+ x=160;
+ z0[i]=1000;
+ if (track->GetProlongation(0,y,z)) z0[i] = z;
+ dca[i] = track->GetD(0,0);
+ //
+ // dca error parrameterezation + pulls
+ //
+ sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->fP4)*(100*track->fP4));
+ if (TMath::Abs(track->fP3)>1) sdcar[i]*=2.5;
+ cdcar[i] = TMath::Exp((TMath::Abs(track->fP4)-0.0106)*525.3);
+ pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
+ pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
+ pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
+ if (track->fTPCr[1]+track->fTPCr[2]+track->fTPCr[3]>0.5) {
+ if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
+ }
+ if (track->fTPCr[4]>0.5) {
+ if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
+ }
+ if (track->fTPCr[0]>0.4) {
+ isPrim[i]=kFALSE; //electron no sigma cut
+ }
+ }
+ //
+ //
+ TStopwatch timer;
+ timer.Start();
+ Int_t ncandidates =0;
+ Int_t nall =0;
+ Int_t ntracks=0;
+ Double_t phase[2][2],radius[2];
+ //
+ // Finf V0s loop
+ //
+ //
+ // //
+ TTreeSRedirector &cstream = *fDebugStreamer;
+ Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
+ AliESDV0MI vertex;
+ Double_t cradius0 = 10*10;
+ Double_t cradius1 = 200*200;
+ Double_t cdist1=3.;
+ Double_t cdist2=4.;
+ Double_t cpointAngle = 0.95;
+ //
+ Double_t delta[2]={10000,10000};
+ for (Int_t i =0;i<nentries;i++){
+ if (sign[i]==0) continue;
+ AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+ if (!track0) continue;
+ cstream<<"Tracks"<<
+ "Tr0.="<<track0<<
+ "dca="<<dca[i]<<
+ "z0="<<z0[i]<<
+ "zvertex="<<zvertex<<
+ "sdcar0="<<sdcar[i]<<
+ "cdcar0="<<cdcar[i]<<
+ "pulldcar0="<<pulldcar[i]<<
+ "pulldcaz0="<<pulldcaz[i]<<
+ "pulldca0="<<pulldca[i]<<
+ "isPrim="<<isPrim[i]<<
+ "\n";
+ //
+ if (track0->fP4<0) continue;
+ if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
+ //
+ if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
+ ntracks++;
+ // debug output
+
+
+ for (Int_t j =0;j<nentries;j++){
+ AliTPCseed * track1 = (AliTPCseed*)array->At(j);
+ if (!track1) continue;
+ if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
+ if (sign[j]*sign[i]>0) continue;
+ if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
+ if (track0->fCircular+track1->fCircular>1) continue; //circling -returning track
+ nall++;
+ //
+ // DCA to prim vertex cut
+ //
+ //
+ delta[0]=10000;
+ delta[1]=10000;
+
+ Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
+ if (npoints<1) continue;
+ Int_t iclosest=0;
+ // cuts on radius
+ if (npoints==1){
+ if (radius[0]<cradius0||radius[0]>cradius1) continue;
+ helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
+ if (delta[0]>cdist1) continue;
+ }
+ else{
+ if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
+ helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
+ helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
+ if (delta[1]<delta[0]) iclosest=1;
+ if (delta[iclosest]>cdist1) continue;
+ }
+ helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
+ if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
+ //
+ Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
+ if (pointAngle<cpointAngle) continue;
+ //
+ Bool_t isGamma = kFALSE;
+ vertex.SetP(*track0); //track0 - plus
+ vertex.SetM(*track1); //track1 - minus
+ vertex.Update(fprimvertex);
+ if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
+ Double_t pointAngle2 = vertex.GetPointAngle();
+ //continue;
+ if (vertex.GetPointAngle()<cpointAngle && (!isGamma)) continue; // point angle cut
+ if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
+ Float_t sigmae = 0.15*0.15;
+ if (vertex.GetRr()<80)
+ sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
+ sigmae+= TMath::Sqrt(sigmae);
+ if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
+ Float_t densb0=0,densb1=0,densa0=0,densa1=0;
+ Int_t row0 = GetRowNumber(vertex.GetRr());
+ if (row0>15){
+ if (vertex.GetDist2()>0.2) continue;
+ densb0 = track0->Density2(0,row0-5);
+ densb1 = track1->Density2(0,row0-5);
+ if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
+ densa0 = track0->Density2(row0+5,row0+40);
+ densa1 = track1->Density2(row0+5,row0+40);
+ if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
+ }
+ else{
+ densa0 = track0->Density2(0,40); //cluster density
+ densa1 = track1->Density2(0,40); //cluster density
+ if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
+ }
+ vertex.SetLab(0,track0->GetLabel());
+ vertex.SetLab(1,track1->GetLabel());
+ vertex.SetChi2After((densa0+densa1)*0.5);
+ vertex.SetChi2Before((densb0+densb1)*0.5);
+ vertex.SetIndex(0,i);
+ vertex.SetIndex(1,j);
+ vertex.SetStatus(1); // TPC v0 candidate
+ vertex.SetRp(track0->fTPCr);
+ vertex.SetRm(track1->fTPCr);
+ tpcv0s->AddLast(new AliESDV0MI(vertex));
+ ncandidates++;
+ {
+ Int_t eventNr = esd->GetEventNumber();
+ Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
+ Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
+ cstream<<"V0"<<
+ "Event="<<eventNr<<
+ "vertex.="<<&vertex<<
+ "Tr0.="<<track0<<
+ "lab0="<<track0->fLab<<
+ "Helix0.="<<&helixes[i]<<
+ "Tr1.="<<track1<<
+ "lab1="<<track1->fLab<<
+ "Helix1.="<<&helixes[j]<<
+ "pointAngle="<<pointAngle<<
+ "pointAngle2="<<pointAngle2<<
+ "dca0="<<dca[i]<<
+ "dca1="<<dca[j]<<
+ "z0="<<z0[i]<<
+ "z1="<<z0[j]<<
+ "zvertex="<<zvertex<<
+ "circular0="<<track0->fCircular<<
+ "circular1="<<track1->fCircular<<
+ "npoints="<<npoints<<
+ "radius0="<<radius[0]<<
+ "delta0="<<delta[0]<<
+ "radius1="<<radius[1]<<
+ "delta1="<<delta[1]<<
+ "radiusm="<<radiusm<<
+ "deltam="<<deltam<<
+ "sdcar0="<<sdcar[i]<<
+ "sdcar1="<<sdcar[j]<<
+ "cdcar0="<<cdcar[i]<<
+ "cdcar1="<<cdcar[j]<<
+ "pulldcar0="<<pulldcar[i]<<
+ "pulldcar1="<<pulldcar[j]<<
+ "pulldcaz0="<<pulldcaz[i]<<
+ "pulldcaz1="<<pulldcaz[j]<<
+ "pulldca0="<<pulldca[i]<<
+ "pulldca1="<<pulldca[j]<<
+ "densb0="<<densb0<<
+ "densb1="<<densb1<<
+ "densa0="<<densa0<<
+ "densa1="<<densa1<<
+ "sigmae="<<sigmae<<
+ "\n";
+ }
+ }
+ }
+ Float_t *quality = new Float_t[ncandidates];
+ Int_t *indexes = new Int_t[ncandidates];
+ Int_t naccepted =0;
+ for (Int_t i=0;i<ncandidates;i++){
+ quality[i] = 0;
+ AliESDV0MI *v0 = (AliESDV0MI*)tpcv0s->At(i);
+ quality[i] = 1./(1.00001-v0->GetPointAngle()); //base point angle
+ // quality[i] /= (0.5+v0->GetDist2());
+ // quality[i] *= v0->GetChi2After(); //density factor
+ Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
+ Int_t index0 = v0->GetIndex(0);
+ Int_t index1 = v0->GetIndex(1);
+ AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
+ if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
+ if (track0->fTPCr[4]>0.9||track1->fTPCr[4]>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
+ }
+
+ TMath::Sort(ncandidates,quality,indexes,kTRUE);
+ //
+ //
+ for (Int_t i=0;i<ncandidates;i++){
+ AliESDV0MI * v0 = (AliESDV0MI*)tpcv0s->At(indexes[i]);
+ if (!v0) continue;
+ Int_t index0 = v0->GetIndex(0);
+ Int_t index1 = v0->GetIndex(1);
+ AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
+ AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
+ if (!track0||!track1) {
+ printf("Bug\n");
+ continue;
+ }
+ Bool_t accept =kTRUE; //default accept
+ Int_t *v0indexes0 = track0->GetV0Indexes();
+ Int_t *v0indexes1 = track1->GetV0Indexes();
+ //
+ Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
+ Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
+ if (v0indexes0[1]!=0) order0 =2;
+ if (v0indexes1[1]!=0) order1 =2;
+ //
+ if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
+ if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
+ //
+ AliESDV0MI * v02 = v0;
+ if (accept){
+ v0->SetOrder(0,order0);
+ v0->SetOrder(1,order1);
+ v0->SetOrder(1,order0+order1);
+ Int_t index = esd->AddV0MI(v0);
+ v02 = esd->GetV0MI(index);
+ v0indexes0[order0]=index;
+ v0indexes1[order1]=index;
+ naccepted++;
+ }
+ {
+ Int_t eventNr = esd->GetEventNumber();
+ cstream<<"V02"<<
+ "Event="<<eventNr<<
+ "vertex.="<<v0<<
+ "vertex2.="<<v02<<
+ "Tr0.="<<track0<<
+ "lab0="<<track0->fLab<<
+ "Tr1.="<<track1<<
+ "lab1="<<track1->fLab<<
+ "dca0="<<dca[index0]<<
+ "dca1="<<dca[index1]<<
+ "order0="<<order0<<
+ "order1="<<order1<<
+ "accept="<<accept<<
+ "quality="<<quality[i]<<
+ "pulldca0="<<pulldca[index0]<<
+ "pulldca1="<<pulldca[index1]<<
+ "index="<<i<<
+ "\n";
+ }
+ }
+
+
+ //
+ //
+ delete []quality;
+ delete []indexes;
+//
+ delete [] isPrim;
+ delete [] pulldca;
+ delete [] pulldcaz;
+ delete [] pulldcar;
+ delete [] cdcar;
+ delete [] sdcar;
+ delete [] dca;
+ //
+ delete[] z0;
+ delete[] alpha;
+ delete[] sign;
+ delete[] helixes;
+ printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
+ timer.Print();
+}
+
Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink)
{
//
//
RemoveUsed2(fSeeds,0.85,0.85,0);
FindKinks(fSeeds,fEvent);
- RemoveUsed2(fSeeds,0.5,0.4,50);
+ RemoveUsed2(fSeeds,0.5,0.4,20);
+ // //
+// // refit short tracks
+// //
+ Int_t nseed=fSeeds->GetEntriesFast();
+// for (Int_t i=0; i<nseed; i++) {
+// AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
+// if (!pt) continue;
+// Int_t nc=t.GetNumberOfClusters();
+// if (nc<15) {
+// delete fSeeds->RemoveAt(i);
+// continue;
+// }
+// if (pt->GetKinkIndexes()[0]!=0) continue; // ignore kink candidates
+// if (nc>100) continue; // hopefully, enough for ITS
+// AliTPCseed *seed2 = new AliTPCseed(*pt);
+// //seed2->Reset(kFALSE);
+// //pt->ResetCovariance();
+// seed2->Modify(1);
+// FollowBackProlongation(*seed2,158);
+// //seed2->Reset(kFALSE);
+// seed2->Modify(10);
+// FollowProlongation(*seed2,0);
+// TTreeSRedirector &cstream = *fDebugStreamer;
+// cstream<<"Crefit"<<
+// "Tr0.="<<pt<<
+// "Tr1.="<<seed2<<
+// "\n";
+// if (seed2->fN>pt->fN){
+// delete fSeeds->RemoveAt(i);
+// fSeeds->AddAt(seed2,i);
+// }else{
+// delete seed2;
+// }
+// }
+// RemoveUsed2(fSeeds,0.6,0.6,50);
+
+// FindV0s(fSeeds,fEvent);
//RemoveDouble(fSeeds,0.2,0.6,11);
- // RemoveUsed(fSeeds,0.5,0.5,6);
//
- Int_t nseed=fSeeds->GetEntriesFast();
Int_t found = 0;
for (Int_t i=0; i<nseed; i++) {
AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
-AliTPCseed::AliTPCseed():AliTPCtrack(){
- //
- fRow=0;
- fRemoval =0;
- for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
- for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
- for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
- fPoints = 0;
- fEPoints = 0;
- fNFoundable =0;
- fNShared =0;
- fRemoval = 0;
- fSort =0;
- fFirstPoint =0;
- fNoCluster =0;
- fBSigned = kFALSE;
- fSeed1 =-1;
- fSeed2 =-1;
- fCurrentCluster =0;
- fCurrentSigmaY2=0;
- fCurrentSigmaZ2=0;
-}
-AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
- //---------------------
- // dummy copy constructor
- //-------------------------
- for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i];
- for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
- fPoints = 0;
- fEPoints = 0;
-}
-AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
- //
- //copy constructor
- fPoints = 0;
- fEPoints = 0;
- fNShared =0;
- // fTrackPoints =0;
- fRemoval =0;
- fSort =0;
- for (Int_t i=0;i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i);
-
- for (Int_t i=0;i<160;i++) {
- fClusterPointer[i] = 0;
- Int_t index = t.GetClusterIndex(i);
- if (index>=-1){
- SetClusterIndex2(i,index);
- }
- else{
- SetClusterIndex2(i,-3);
- }
- }
- fFirstPoint =0;
- fNoCluster =0;
- fBSigned = kFALSE;
- fSeed1 =-1;
- fSeed2 =-1;
- fCurrentCluster =0;
- fCurrentSigmaY2=0;
- fCurrentSigmaZ2=0;
-}
-/*
-AliTPCseed::AliTPCseed(const AliTPCtrack &t, Double_t a):AliTPCtrack(t,a){
- //
- //copy constructor
- fRow=0;
- for (Int_t i=0;i<160;i++) {
- fClusterPointer[i] = 0;
- Int_t index = t.GetClusterIndex(i);
- SetClusterIndex2(i,index);
- }
- for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
- fPoints = 0;
- fEPoints = 0;
- fNFoundable =0;
- fNShared =0;
- // fTrackPoints =0;
- fRemoval =0;
- fSort = 0;
- fFirstPoint =0;
- fNoCluster =0;
- fBSigned = kFALSE;
- fSeed1 =-1;
- fSeed2 =-1;
- fCurrentCluster =0;
- fCurrentSigmaY2=0;
- fCurrentSigmaZ2=0;
-
-}
-*/
-
-AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
- Double_t xr, Double_t alpha):
- AliTPCtrack(index, xx, cc, xr, alpha) {
- //
- //
- //constructor
- fRow =0;
- for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
- for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
- for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
- fPoints = 0;
- fEPoints = 0;
- fNFoundable =0;
- fNShared = 0;
- // fTrackPoints =0;
- fRemoval =0;
- fSort =0;
- fFirstPoint =0;
- // fHelixIn = new TClonesArray("AliHelix",0);
- //fHelixOut = new TClonesArray("AliHelix",0);
- fNoCluster =0;
- fBSigned = kFALSE;
- fSeed1 =-1;
- fSeed2 =-1;
- fCurrentCluster =0;
- fCurrentSigmaY2=0;
- fCurrentSigmaZ2=0;
-}
-
-AliTPCseed::~AliTPCseed(){
- //
- // destructor
- if (fPoints) delete fPoints;
- fPoints =0;
- if (fEPoints) delete fEPoints;
- fEPoints = 0;
- fNoCluster =0;
-}
-
-AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
-{
- //
- //
- return &fTrackPoints[i];
-}
-
-void AliTPCseed::RebuildSeed()
-{
- //
- // rebuild seed to be ready for storing
- AliTPCclusterMI cldummy;
- cldummy.SetQ(0);
- AliTPCTrackPoint pdummy;
- pdummy.GetTPoint().fIsShared = 10;
- for (Int_t i=0;i<160;i++){
- AliTPCclusterMI * cl0 = fClusterPointer[i];
- AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
- if (cl0){
- trpoint->GetTPoint() = *(GetTrackPoint(i));
- trpoint->GetCPoint() = *cl0;
- trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
- }
- else{
- *trpoint = pdummy;
- trpoint->GetCPoint()= cldummy;
- }
-
- }
-
-}
-
-
-Double_t AliTPCseed::GetDensityFirst(Int_t n)
-{
- //
- //
- // return cluster for n rows bellow first point
- Int_t nfoundable = 1;
- Int_t nfound = 1;
- for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
- Int_t index = GetClusterIndex2(i);
- if (index!=-1) nfoundable++;
- if (index>0) nfound++;
- }
- if (nfoundable<n) return 0;
- return Double_t(nfound)/Double_t(nfoundable);
-
-}
-
-
-void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
-{
- // get cluster stat. on given region
- //
- found = 0;
- foundable = 0;
- shared =0;
- for (Int_t i=first;i<last; i++){
- Int_t index = GetClusterIndex2(i);
- if (index!=-1) foundable++;
- if (fClusterPointer[i]) {
- found++;
- }
- else
- continue;
-
- if (fClusterPointer[i]->IsUsed(10)) {
- shared++;
- continue;
- }
- if (!plus2) continue; //take also neighborhoud
- //
- if ( (i>0) && fClusterPointer[i-1]){
- if (fClusterPointer[i-1]->IsUsed(10)) {
- shared++;
- continue;
- }
- }
- if ( fClusterPointer[i+1]){
- if (fClusterPointer[i+1]->IsUsed(10)) {
- shared++;
- continue;
- }
- }
-
- }
- //if (shared>found){
- //Error("AliTPCseed::GetClusterStatistic","problem\n");
- //}
-}
-
-//_____________________________________________________________________________
-void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
- //-----------------------------------------------------------------
- // This funtion calculates dE/dX within the "low" and "up" cuts.
- //-----------------------------------------------------------------
-
- Float_t amp[200];
- Float_t angular[200];
- Float_t weight[200];
- Int_t index[200];
- //Int_t nc = 0;
- // TClonesArray & arr = *fPoints;
- Float_t meanlog = 100.;
-
- Float_t mean[4] = {0,0,0,0};
- Float_t sigma[4] = {1000,1000,1000,1000};
- Int_t nc[4] = {0,0,0,0};
- Float_t norm[4] = {1000,1000,1000,1000};
- //
- //
- fNShared =0;
-
- for (Int_t of =0; of<4; of++){
- for (Int_t i=of+i1;i<i2;i+=4)
- {
- Int_t index = fIndex[i];
- if (index<0||index&0x8000) continue;
-
- //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
- AliTPCTrackerPoint * point = GetTrackPoint(i);
- //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
- //AliTPCTrackerPoint * pointp = 0;
- //if (i<159) pointp = GetTrackPoint(i+1);
-
- if (point==0) continue;
- AliTPCclusterMI * cl = fClusterPointer[i];
- if (cl==0) continue;
- if (onlyused && (!cl->IsUsed(10))) continue;
- if (cl->IsUsed(11)) {
- fNShared++;
- continue;
- }
- Int_t type = cl->GetType();
- //if (point->fIsShared){
- // fNShared++;
- // continue;
- //}
- //if (pointm)
- // if (pointm->fIsShared) continue;
- //if (pointp)
- // if (pointp->fIsShared) continue;
-
- if (type<0) continue;
- //if (type>10) continue;
- //if (point->GetErrY()==0) continue;
- //if (point->GetErrZ()==0) continue;
-
- //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
- //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
- //if ((ddy*ddy+ddz*ddz)>10) continue;
-
-
- // if (point->GetCPoint().GetMax()<5) continue;
- if (cl->GetMax()<5) continue;
- Float_t angley = point->GetAngleY();
- Float_t anglez = point->GetAngleZ();
-
- Float_t rsigmay2 = point->GetSigmaY();
- Float_t rsigmaz2 = point->GetSigmaZ();
- /*
- Float_t ns = 1.;
- if (pointm){
- rsigmay += pointm->GetTPoint().GetSigmaY();
- rsigmaz += pointm->GetTPoint().GetSigmaZ();
- ns+=1.;
- }
- if (pointp){
- rsigmay += pointp->GetTPoint().GetSigmaY();
- rsigmaz += pointp->GetTPoint().GetSigmaZ();
- ns+=1.;
- }
- rsigmay/=ns;
- rsigmaz/=ns;
- */
-
- Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
-
- Float_t ampc = 0; // normalization to the number of electrons
- if (i>64){
- // ampc = 1.*point->GetCPoint().GetMax();
- ampc = 1.*cl->GetMax();
- //ampc = 1.*point->GetCPoint().GetQ();
- // AliTPCClusterPoint & p = point->GetCPoint();
- // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
- // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
- //Float_t dz =
- // TMath::Abs( Int_t(iz) - iz + 0.5);
- //ampc *= 1.15*(1-0.3*dy);
- //ampc *= 1.15*(1-0.3*dz);
- // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
- //ampc *=zfactor;
- }
- else{
- //ampc = 1.0*point->GetCPoint().GetMax();
- ampc = 1.0*cl->GetMax();
- //ampc = 1.0*point->GetCPoint().GetQ();
- //AliTPCClusterPoint & p = point->GetCPoint();
- // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
- //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
- //Float_t dz =
- // TMath::Abs( Int_t(iz) - iz + 0.5);
-
- //ampc *= 1.15*(1-0.3*dy);
- //ampc *= 1.15*(1-0.3*dz);
- // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
- //ampc *=zfactor;
-
- }
- ampc *= 2.0; // put mean value to channel 50
- //ampc *= 0.58; // put mean value to channel 50
- Float_t w = 1.;
- // if (type>0) w = 1./(type/2.-0.5);
- // Float_t z = TMath::Abs(cl->GetZ());
- if (i<64) {
- ampc /= 0.6;
- //ampc /= (1+0.0008*z);
- } else
- if (i>128){
- ampc /=1.5;
- //ampc /= (1+0.0008*z);
- }else{
- //ampc /= (1+0.0008*z);
- }
-
- if (type<0) { //amp at the border - lower weight
- // w*= 2.;
-
- continue;
- }
- if (rsigma>1.5) ampc/=1.3; // if big backround
- amp[nc[of]] = ampc;
- angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
- weight[nc[of]] = w;
- nc[of]++;
- }
-
- TMath::Sort(nc[of],amp,index,kFALSE);
- Float_t sumamp=0;
- Float_t sumamp2=0;
- Float_t sumw=0;
- //meanlog = amp[index[Int_t(nc[of]*0.33)]];
- meanlog = 50;
- for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
- Float_t ampl = amp[index[i]]/angular[index[i]];
- ampl = meanlog*TMath::Log(1.+ampl/meanlog);
- //
- sumw += weight[index[i]];
- sumamp += weight[index[i]]*ampl;
- sumamp2 += weight[index[i]]*ampl*ampl;
- norm[of] += angular[index[i]]*weight[index[i]];
- }
- if (sumw<1){
- SetdEdx(0);
- }
- else {
- norm[of] /= sumw;
- mean[of] = sumamp/sumw;
- sigma[of] = sumamp2/sumw-mean[of]*mean[of];
- if (sigma[of]>0.1)
- sigma[of] = TMath::Sqrt(sigma[of]);
- else
- sigma[of] = 1000;
-
- mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
- //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
- //mean *=(1-0.1*(norm-1.));
- }
- }
-
- Float_t dedx =0;
- fSdEdx =0;
- fMAngular =0;
- // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
- // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
-
-
- // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
- // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
-
- Int_t norm2 = 0;
- Int_t norm3 = 0;
- for (Int_t i =0;i<4;i++){
- if (nc[i]>2&&nc[i]<1000){
- dedx += mean[i] *nc[i];
- fSdEdx += sigma[i]*(nc[i]-2);
- fMAngular += norm[i] *nc[i];
- norm2 += nc[i];
- norm3 += nc[i]-2;
- }
- fDEDX[i] = mean[i];
- fSDEDX[i] = sigma[i];
- fNCDEDX[i]= nc[i];
- }
-
- if (norm3>0){
- dedx /=norm2;
- fSdEdx /=norm3;
- fMAngular/=norm2;
- }
- else{
- SetdEdx(0);
- return;
- }
- // Float_t dedx1 =dedx;
- /*
- dedx =0;
- for (Int_t i =0;i<4;i++){
- if (nc[i]>2&&nc[i]<1000){
- mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
- dedx += mean[i] *nc[i];
- }
- fDEDX[i] = mean[i];
- }
- dedx /= norm2;
- */
-
-
- SetdEdx(dedx);
-
- //mi deDX
-
-
-
- //Very rough PID
- Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
-
- if (p<0.6) {
- if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
- if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
- SetMass(0.93827); return;
- }
-
- if (p<1.2) {
- if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
- SetMass(0.93827); return;
- }
-
- SetMass(0.13957); return;
-
-}
-
-
-
-/*
-
-
-void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
- //-----------------------------------------------------------------
- // This funtion calculates dE/dX within the "low" and "up" cuts.
- //-----------------------------------------------------------------
-
- Float_t amp[200];
- Float_t angular[200];
- Float_t weight[200];
- Int_t index[200];
- Bool_t inlimit[200];
- for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
- for (Int_t i=0;i<200;i++) amp[i]=10000;
- for (Int_t i=0;i<200;i++) angular[i]= 1;;
-
-
- //
- Float_t meanlog = 100.;
- Int_t indexde[4]={0,64,128,160};
+// AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
+// {
+// //
+// //
+// return &fTrackPoints[i];
+// }
- Float_t amean =0;
- Float_t asigma =0;
- Float_t anc =0;
- Float_t anorm =0;
-
- Float_t mean[4] = {0,0,0,0};
- Float_t sigma[4] = {1000,1000,1000,1000};
- Int_t nc[4] = {0,0,0,0};
- Float_t norm[4] = {1000,1000,1000,1000};
- //
- //
- fNShared =0;
-
- // for (Int_t of =0; of<3; of++){
- // for (Int_t i=indexde[of];i<indexde[of+1];i++)
- for (Int_t i =0; i<160;i++)
- {
- AliTPCTrackPoint * point = GetTrackPoint(i);
- if (point==0) continue;
- if (point->fIsShared){
- fNShared++;
- continue;
- }
- Int_t type = point->GetCPoint().GetType();
- if (type<0) continue;
- if (point->GetCPoint().GetMax()<5) continue;
- Float_t angley = point->GetTPoint().GetAngleY();
- Float_t anglez = point->GetTPoint().GetAngleZ();
- Float_t rsigmay = point->GetCPoint().GetSigmaY();
- Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
- Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
-
- Float_t ampc = 0; // normalization to the number of electrons
- if (i>64){
- ampc = point->GetCPoint().GetMax();
- }
- else{
- ampc = point->GetCPoint().GetMax();
- }
- ampc *= 2.0; // put mean value to channel 50
- // ampc *= 0.565; // put mean value to channel 50
-
- Float_t w = 1.;
- Float_t z = TMath::Abs(point->GetCPoint().GetZ());
- if (i<64) {
- ampc /= 0.63;
- } else
- if (i>128){
- ampc /=1.51;
- }
- if (type<0) { //amp at the border - lower weight
- continue;
- }
- if (rsigma>1.5) ampc/=1.3; // if big backround
- angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
- amp[i] = ampc/angular[i];
- weight[i] = w;
- anc++;
- }
- TMath::Sort(159,amp,index,kFALSE);
- for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
- inlimit[index[i]] = kTRUE; // take all clusters
- }
-
- // meanlog = amp[index[Int_t(anc*0.3)]];
- meanlog =10000.;
- for (Int_t of =0; of<3; of++){
- Float_t sumamp=0;
- Float_t sumamp2=0;
- Float_t sumw=0;
- for (Int_t i=indexde[of];i<indexde[of+1];i++)
- {
- if (inlimit[i]==kFALSE) continue;
- Float_t ampl = amp[i];
- ///angular[i];
- ampl = meanlog*TMath::Log(1.+ampl/meanlog);
- //
- sumw += weight[i];
- sumamp += weight[i]*ampl;
- sumamp2 += weight[i]*ampl*ampl;
- norm[of] += angular[i]*weight[i];
- nc[of]++;
- }
- if (sumw<1){
- SetdEdx(0);
- }
- else {
- norm[of] /= sumw;
- mean[of] = sumamp/sumw;
- sigma[of] = sumamp2/sumw-mean[of]*mean[of];
- if (sigma[of]>0.1)
- sigma[of] = TMath::Sqrt(sigma[of]);
- else
- sigma[of] = 1000;
- mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
- }
- }
-
- Float_t dedx =0;
- fSdEdx =0;
- fMAngular =0;
- //
- Int_t norm2 = 0;
- Int_t norm3 = 0;
- Float_t www[3] = {12.,14.,17.};
- //Float_t www[3] = {1.,1.,1.};
-
- for (Int_t i =0;i<3;i++){
- if (nc[i]>2&&nc[i]<1000){
- dedx += mean[i] *nc[i]*www[i]/sigma[i];
- fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
- fMAngular += norm[i] *nc[i];
- norm2 += nc[i]*www[i]/sigma[i];
- norm3 += (nc[i]-2)*www[i]/sigma[i];
- }
- fDEDX[i] = mean[i];
- fSDEDX[i] = sigma[i];
- fNCDEDX[i]= nc[i];
- }
-
- if (norm3>0){
- dedx /=norm2;
- fSdEdx /=norm3;
- fMAngular/=norm2;
- }
- else{
- SetdEdx(0);
- return;
- }
- // Float_t dedx1 =dedx;
-
- dedx =0;
- Float_t norm4 = 0;
- for (Int_t i =0;i<3;i++){
- if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
- //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
- dedx += mean[i] *(nc[i])/(sigma[i]);
- norm4 += (nc[i])/(sigma[i]);
- }
- fDEDX[i] = mean[i];
- }
- if (norm4>0) dedx /= norm4;
-
-
-
- SetdEdx(dedx);
-
- //mi deDX
-
-}
-
-*/
// Origin:
//-------------------------------------------------------
+#include <TError.h>
#include "AliTracker.h"
-#include "AliTPCtrack.h"
-#include "AliComplexCluster.h"
+#include "AliTPCreco.h"
+#include "AliPID.h"
+
class TFile;
class AliTPCParam;
class AliESD;
class TTree;
class AliESDkink;
-
-class AliTPCseed : public AliTPCtrack {
- friend class AliTPCtrackerMI;
- public:
- AliTPCseed();
- virtual ~AliTPCseed();
- AliTPCseed(const AliTPCtrack &t);
- AliTPCseed(const AliTPCseed &s);
- //AliTPCseed(const AliTPCseed &t, Double_t a);
- AliTPCseed(UInt_t index, const Double_t xx[5],
- const Double_t cc[15], Double_t xr, Double_t alpha);
- Int_t Compare(const TObject *o) const;
- void Reset(Bool_t all = kTRUE);
- Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const;
- virtual Double_t GetPredictedChi2(const AliTPCclusterMI *cluster) const;
- virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i);
- AliTPCTrackerPoint * GetTrackPoint(Int_t i);
- void RebuildSeed(); // rebuild seed to be ready for storing
- Double_t GetDensityFirst(Int_t n);
- Double_t GetSigma2C() const {return fC44;};
- void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2);
-
- void Modify(Double_t factor);
- void SetClusterIndex2(Int_t row, Int_t index) {
- fIndex[row] = index;
- }
- Int_t GetClusterIndex2(Int_t row) const {
- return fIndex[row];
- }
- Int_t GetClusterSector(Int_t row) const {
- Int_t pica = -1;
- if (fIndex[row]>=0) pica = ((fIndex[row]&0xff000000)>>24);
- return pica;
- }
-
- void SetErrorY2(Float_t sy2){fErrorY2=sy2;}
- void SetErrorZ2(Float_t sz2){fErrorZ2=sz2;}
- void CookdEdx(Double_t low=0.05, Double_t up=0.70, Int_t i1=0, Int_t i2=159, Bool_t onlyused = kFALSE);
- // void CookdEdx2(Double_t low=0.05, Double_t up=0.70);
- Bool_t IsActive() const { return !(fRemoval);}
- void Desactivate(Int_t reason){ fRemoval = reason;}
- //
- //
- private:
- AliTPCseed & operator = (const AliTPCseed &);
- AliESDtrack * fEsd; //!
- AliTPCclusterMI* fClusterPointer[160]; //! array of cluster pointers -
- TClonesArray * fPoints; // array with points along the track
- TClonesArray * fEPoints; // array with exact points - calculated in special macro not used in tracking
- //---CURRENT VALUES
- Int_t fRow; //!current row number
- Int_t fSector; //!current sector number
- Int_t fRelativeSector; //! index of current relative sector
- Float_t fCurrentSigmaY2; //!expected current cluster sigma Y
- Float_t fCurrentSigmaZ2; //!expected current cluster sigma Z
- Float_t fErrorY2; //!sigma of current cluster
- Float_t fErrorZ2; //!sigma of current cluster
- AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation
- Int_t fCurrentClusterIndex1; //! index of the current cluster
- Bool_t fInDead; //! indicate if the track is in dead zone
- Bool_t fIsSeeding; //!indicates if it is proces of seeading
- Int_t fNoCluster; //!indicates number of rows without clusters
- Int_t fSort; //!indicate criteria for sorting
- Bool_t fBSigned; //indicates that clusters of this trackes are signed to be used
- //
- //
- Float_t fDEDX[4]; // dedx according padrows
- Float_t fSDEDX[4]; // sdedx according padrows
- Int_t fNCDEDX[4]; // number of clusters for dedx measurment
- //
- Int_t fSeedType; //seeding type
- Int_t fSeed1; //first row for seeding
- Int_t fSeed2; //last row for seeding
- Int_t fOverlapLabels[12]; //track labels and the length of the overlap
- Float_t fMAngular; // mean angular factor
- AliTPCTrackerPoint fTrackPoints[160]; //!track points - array track points
-
- ClassDef(AliTPCseed,1)
-};
-
-
+class TTreeSRedirector;
class AliTPCtrackerMI : public AliTracker {
void DeleteSeeds();
void SetDebug(Int_t debug){ fDebug = debug;}
void FindKinks(TObjArray * array, AliESD * esd);
+ void FindV0s(TObjArray * array, AliESD * esd);
void UpdateKinkQualityM(AliTPCseed * seed);
void UpdateKinkQualityD(AliTPCseed * seed);
Int_t CheckKinkPoint(AliTPCseed*seed, AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink);
Double_t fYMax[200]; // max y for given pad row
Double_t fPadLength[200]; // max y for given pad row
const AliTPCParam *fParam; //pointer to the parameters
+ TTreeSRedirector *fDebugStreamer; //!debug streamer
ClassDef(AliTPCtrackerMI,1)
};
AliClustersArray.cxx AliTPCClustersArray.cxx \
AliTPCclusterer.cxx AliTPCclustererMI.cxx \
AliTPCtrack.cxx AliTPCtracker.cxx \
- AliTPCpolyTrack.cxx AliTPCtrackerMI.cxx \
+ AliTPCpolyTrack.cxx AliTPCseed.cxx AliTPCtrackerMI.cxx \
AliTPCPid.cxx AliTPCtrackPid.cxx AliTPCpidESD.cxx \
AliTPCReconstructor.cxx
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff