track labels added
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITStrackerUpgrade.cxx
index 86aee2d7a469b4eaf922c07857a66ed2347d832c..c04ac12728583c0d335c4196349ad79b9097651a 100644 (file)
@@ -71,6 +71,7 @@ ClassImp(AliITStrackerUpgrade)
                                               fMinNPoints(0),
                                               fMinQ(0.),
                                               fLayers(0),
+                                               fSegmentation(0x0),
                                               fCluLayer(0),
                                               fCluCoord(0){
     // Default constructor
@@ -101,6 +102,7 @@ AliITStrackerUpgrade::AliITStrackerUpgrade(const AliITStrackerUpgrade& tracker):
                                                                                fMinNPoints(tracker.fMinNPoints),
                                                                                fMinQ(tracker.fMinQ),
                                                                                fLayers(tracker.fLayers),
+                                                                               fSegmentation(tracker.fSegmentation),
                                                                                fCluLayer(tracker.fCluLayer),
                                                                                fCluCoord(tracker.fCluCoord) {
   // Copy constructor
@@ -166,6 +168,7 @@ AliITStrackerUpgrade::~AliITStrackerUpgrade(){
     delete [] fCluCoord;
   }
   
+ if(fSegmentation) delete fSegmentation;
 }
 
 //____________________________________________________________________________
@@ -224,37 +227,37 @@ void AliITStrackerUpgrade::Init(){
   fCluLayer = 0;
   fCluCoord = 0;
   fMinNPoints = 3;
-  AliITSsegmentationUpgrade *segmentation2 = 0x0;
-  if(!segmentation2)
-    segmentation2 = new AliITSsegmentationUpgrade();
   for(Int_t layer=0; layer<6; layer++){
     Double_t p=0.;
     Double_t zC= 0.;
     fLayers[layer] = new AliITSlayerUpgrade(p,zC);
   }
 
+  fSegmentation = new AliITSsegmentationUpgrade();
        
 }
 //_______________________________________________________________________
 Int_t AliITStrackerUpgrade::LoadClusters(TTree *clusTree){ 
+//
+// Load clusters for tracking
+// 
+
   TClonesArray statITSCluster("AliITSRecPoint");
-  for(Int_t layer=0; layer<6; layer++){
-    TClonesArray *ITSCluster = &statITSCluster;
+  TClonesArray *ITSCluster = &statITSCluster;
+
     TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
     if (!itsClusterBranch){
-      printf("can't get the branch with the ITS clusters ! \n");
+      AliError("can't get the branch with the ITS clusters ! \n");
       return 1;
     }
     itsClusterBranch->SetAddress(&ITSCluster);
-    if (!clusTree->GetEvent(layer))    continue;
+    clusTree->GetEvent(0);
     Int_t nCluster = ITSCluster->GetEntriesFast();
     for(Int_t i=0; i<nCluster; i++){
       AliITSRecPoint *recp = (AliITSRecPoint*)ITSCluster->UncheckedAt(i);
-
-      fLayers[layer]->InsertCluster(new AliITSRecPoint(*recp));
+      fLayers[recp->GetLayer()]->InsertCluster(new AliITSRecPoint(*recp));
     }//loop clusters
-  }//loop layer
-   
+
   SetClusterTree(clusTree);
   return 0;
 }
@@ -357,6 +360,7 @@ Int_t AliITStrackerUpgrade::FindTracks(AliESDEvent* event,Bool_t useAllClusters)
     TClonesArray &clulay = *fCluLayer[ilay];
     TClonesArray &clucoo = *fCluCoord[ilay];
     if (!ForceSkippingOfLayer(ilay)){
+    AliDebug(2,Form("number of clusters in layer %i : %i",ilay,fLayers[ilay]->GetNumberOfClusters()));
       for(Int_t cli=0;cli<fLayers[ilay]->GetNumberOfClusters();cli++){
         AliITSRecPoint* cls = (AliITSRecPoint*)fLayers[ilay]->GetCluster(cli);
         if(cls->TestBit(kSAflag)==kTRUE) continue;
@@ -370,7 +374,7 @@ Int_t AliITStrackerUpgrade::FindTracks(AliESDEvent* event,Bool_t useAllClusters)
         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
         dmar[ilay]++;
       }
-    }
+    } else AliDebug(2,Form("Force skipping layer %i",ilay));
   }
 
 
@@ -615,10 +619,7 @@ AliITStrackV2* AliITStrackerUpgrade::FitTrack(AliITStrackSA* tr,Double_t *primar
               Double_t yclu1 = p1->GetY();
               Double_t zclu1 = p1->GetZ();
              layer=p1->GetLayer();
-             AliITSsegmentationUpgrade *segmentation = 0x0;
-             if(!segmentation)
-                segmentation = new AliITSsegmentationUpgrade();
-             radius=segmentation->GetRadius(layer);
+             fSegmentation->GetRadius(layer);
              Double_t cv=0.,tgl2=0.,phi2=0.;
              Int_t cln1=mrk1;
              AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
@@ -647,7 +648,7 @@ AliITStrackV2* AliITStrackerUpgrade::FitTrack(AliITStrackSA* tr,Double_t *primar
              xz[0]= p1->GetDetLocalX(); 
              xz[1]= p1->GetDetLocalZ();   
              Bool_t check2;
-             check2 = segmentation->DetToGlobal(layer,xz[0], xz[1],x,y,z);
+             check2 = fSegmentation->DetToGlobal(layer,xz[0], xz[1],x,y,z);
 
 
               Double_t phiclrad, phicldeg;
@@ -696,14 +697,14 @@ AliITStrackV2* AliITStrackerUpgrade::FitTrack(AliITStrackSA* tr,Double_t *primar
              // Propagate inside the innermost layer with a cluster 
              if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
                Double_t rt=0.;
-               rt=segmentation->GetRadius(5);
+               rt=fSegmentation->GetRadius(5);
                if(RefitAtBase(rt,&ot,inx)){ //fit from layer 1 to layer 6
                  AliITStrackMI otrack2(ot);
                  otrack2.ResetCovariance(10.); 
                  otrack2.ResetClusters();
                  
                  //fit from layer 6 to layer 1
-                 if(RefitAtBase(/*AliITSRecoParam::GetrInsideSPD1()*/segmentation->GetRadius(0),&otrack2,inx)) {//check clind
+                 if(RefitAtBase(/*AliITSRecoParam::GetrInsideSPD1()*/fSegmentation->GetRadius(0),&otrack2,inx)) {//check clind
                    new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
                    new(arrSA[nFoundTracks]) AliITStrackSA(trac);
                    ++nFoundTracks;
@@ -730,6 +731,12 @@ AliITStrackV2* AliITStrackerUpgrade::FitTrack(AliITStrackSA* tr,Double_t *primar
   if(otrack==0) {
     return 0;
   }
+
+  CookLabel(otrack,0.); //MI change - to see fake ratio
+  Int_t label=FindLabel(otrack);
+  otrack->SetLabel(label);
+
+
   Int_t indexc[6];
   for(Int_t i=0;i<6;i++) indexc[i]=0;
   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
@@ -772,14 +779,14 @@ void AliITStrackerUpgrade::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_
 Int_t AliITStrackerUpgrade::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
   //function used to to find the clusters associated to the track
 
-  AliITSsegmentationUpgrade *segmentation = 0x0;
-  if(!segmentation)
-    segmentation = new AliITSsegmentationUpgrade();
-
-  if(ForceSkippingOfLayer(layer)) return 0;
+  AliDebug(2,"Starting...");
+  if(ForceSkippingOfLayer(layer)) {
+  AliDebug(2,Form("Forcing skipping of layer %i. Exiting",layer));
+  return 0;
+  }
 
   Int_t nc=0;
-  Double_t r=segmentation->GetRadius(layer);
+  Double_t r=fSegmentation->GetRadius(layer);
   if(pflag==1){      
     Float_t cx1,cx2,cy1,cy2;
     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
@@ -792,6 +799,7 @@ Int_t AliITStrackerUpgrade::SearchClusters(Int_t layer,Double_t phiwindow,Double
 
  
   Int_t ncl = fCluLayer[layer]->GetEntries();
+  AliDebug(2,Form(" Number of clusters %i in layer %i.",ncl,layer));
   for (Int_t index=0; index<ncl; index++) {
     AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
     if (!c) continue;
@@ -955,49 +963,61 @@ Int_t AliITStrackerUpgrade::FindTrackLowChiSquare() const {
 }
 
 //__________________________________________________________
-/*Int_t AliITStrackerUpgrade::FindLabel(Int_t l0, Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5){
+Int_t AliITStrackerUpgrade::FindLabel(AliITStrackV2* track) const {
+  //
 
-//function used to determine the track label
-  
-Int_t lb[6] = {l0,l1,l2,l3,l4,l5};
-Int_t aa[6]={1,1,1,1,1,1};
-Int_t ff=0; 
-Int_t ll=0;
-Int_t k=0;Int_t w=0;Int_t num=6;
-for(Int_t i=5;i>=0;i--) if(lb[i]==-1) num=i;
-  
-while(k<num){
-  
-for(Int_t i=k+1;i<num;i++){
-    
-if(lb[k]==lb[i] && aa[k]!=0){
-      
-aa[k]+=1;
-aa[i]=0;
-}
-}
-k++;
-}
+  Int_t labl[AliITSgeomTGeo::kNLayers][3];
+  Int_t cnts[AliITSgeomTGeo::kNLayers][3];
+  for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
+    for(Int_t k=0;k<3;k++){
+      labl[j][k]=-2;
+      cnts[j][k]=1;
+    }
+  }
+  Int_t iNotLabel=0;
+  for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
+    Int_t indexc = track->GetClusterIndex(i);
+    AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
+    AliDebug(2,Form(" cluster index %i  ",indexc));
+    Int_t iLayer=cl->GetLayer();
+    for(Int_t k=0;k<3;k++){
+      labl[iLayer][k]=cl->GetLabel(k);
+      if(labl[iLayer][k]<0) iNotLabel++;
+    }
+  }
+  if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
+
+  for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
+    for(Int_t j2=0; j2<j1;  j2++){
+      for(Int_t k1=0; k1<3; k1++){
+        for(Int_t k2=0; k2<3; k2++){
+          if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
+            cnts[j2][k2]++;
+            cnts[j1][k1]=0;
+          }
+        }
+      }
+    }
+  }
 
-while(w<num){
-for(Int_t j=0;j<6;j++){
-if(aa[w]<aa[j]){
-ff=aa[w];
-aa[w]=aa[j];
-aa[j]=ff;
-ll=lb[w];
-lb[w]=lb[j];
-lb[j]=ll;
-}
-}
-w++;
-}
-  
-if(num<1) return -1; 
-return lb[num-1];
+
+  Int_t cntMax=0;
+  Int_t label=-1;
+  for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
+    for(Int_t k=0;k<3;k++){
+      if(cnts[j][k]>cntMax && labl[j][k]>=0){
+        cntMax=cnts[j][k];
+        label=labl[j][k];
+      }
+    }
+  }
+
+  Int_t lflag=0;
+  for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
+  if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
+  if(lflag<track->GetNumberOfClusters()) label = -label;
+  return label;
 }
-*/
 //_____________________________________________________________________________
 /*Int_t AliITStrackerUpgrade::Label(Int_t gl1, Int_t gl2, Int_t gl3, Int_t gl4, Int_t gl5, Int_t gl6,Int_t gl7, Int_t gl8, Int_t gl9, Int_t gl10,Int_t gl11,
   Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t minNPoints){
@@ -1105,13 +1125,10 @@ void AliITStrackerUpgrade::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double
   Double_t xz[2];
   xz[0]= cl->GetDetLocalX(); 
   xz[1]= cl->GetDetLocalZ() ; 
-  AliITSsegmentationUpgrade *segmentation2 = 0x0;
-  if(!segmentation2)
-    segmentation2 = new AliITSsegmentationUpgrade();
   Bool_t check2;
   Int_t ilayer;
   ilayer = cl->GetLayer();
-  check2 = segmentation2->DetToGlobal(ilayer,xz[0], xz[1],x,y,z);
+  check2 = fSegmentation->DetToGlobal(ilayer,xz[0], xz[1],x,y,z);
 
   if(x!=0 && y!=0)  
     phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
@@ -1173,9 +1190,6 @@ Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackMI *track,
   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
     index[k]=clusters[k];
   }
-  AliITSsegmentationUpgrade *segmentation = 0x0;
-  if(!segmentation)
-    segmentation = new AliITSsegmentationUpgrade();
 
   ULong_t trStatus=0;
   if(track->GetESDtrack()) trStatus=track->GetStatus();
@@ -1184,7 +1198,7 @@ Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackMI *track,
     innermostlayer=5;
     Double_t drphi = TMath::Abs(track->GetD(0.,0.));
     for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
-      if( (drphi < (segmentation->GetRadius(innermostlayer)+1.)) ||
+      if( (drphi < (fSegmentation->GetRadius(innermostlayer)+1.)) ||
           index[innermostlayer] >= 0 ) break;
     }
     AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
@@ -1202,7 +1216,7 @@ Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackMI *track,
 
   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
     Double_t r=0.;
-    r=segmentation->GetRadius(ilayer);
+    r=fSegmentation->GetRadius(ilayer);
 
     if (step<0 && xx>r){
       break;
@@ -1244,7 +1258,6 @@ Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackMI *track,
     Double_t alpha= (ladder*18.+9.)*TMath::Pi()/180.;
     
     if (!track->Propagate(alpha,r)) {
       return kFALSE;
     }
 
@@ -1263,7 +1276,7 @@ Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackMI *track,
          clAcc=cl;                                                              
          maxchi2=chi2;                                                          
        } else {                                                                 
-         return kFALSE;                                                         
+          return kFALSE;                                                         
        }                                                                        
       }                                                                          
     }              
@@ -1287,7 +1300,7 @@ Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackMI *track,
 
   if (!track->PropagateTo(xx,0.,0.)){
     return kFALSE;
-  }
+  } 
   return kTRUE;
 }
 
@@ -1371,10 +1384,7 @@ Int_t AliITStrackerUpgrade::CorrectForLayerMaterial(AliITStrackMI *t,
   }
   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
   Double_t r = 0.;
-  AliITSsegmentationUpgrade *segmentation = 0x0;
-  if(!segmentation)
-    segmentation = new AliITSsegmentationUpgrade();
-  r=segmentation->GetRadius(layerindex);
+  r=fSegmentation->GetRadius(layerindex);
   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
   Double_t xToGo;
@@ -1395,7 +1405,7 @@ Int_t AliITStrackerUpgrade::CorrectForLayerMaterial(AliITStrackMI *t,
   switch(mode) {
   case 1:
     x0=21.82;
-    xOverX0 = segmentation->GetThickness(layerindex)/x0;
+    xOverX0 = fSegmentation->GetThickness(layerindex)/x0;
 
     lengthTimesMeanDensity = xOverX0*x0;
     lengthTimesMeanDensity *= dir;
@@ -1560,6 +1570,7 @@ Int_t AliITStrackerUpgrade::PropagateBack(AliESDEvent *event) {
     if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
     if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 
+/*
     AliITStrackMI *t=0;
     try {
       t=new AliITStrackMI(*esd);
@@ -1568,6 +1579,8 @@ Int_t AliITStrackerUpgrade::PropagateBack(AliESDEvent *event) {
       delete t;
       continue;
     }
+*/
+    AliITStrackMI *t = new AliITStrackMI(*esd);
     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 
     ResetTrackToFollow(*t);
@@ -1609,7 +1622,7 @@ Int_t AliITStrackerUpgrade::PropagateBack(AliESDEvent *event) {
       //fTrackToFollow.CookdEdx();
       //CookLabel(&fTrackToFollow,0.); //For comparison only
       fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
-      UseClusters(&fTrackToFollow);
+      //UseClusters(&fTrackToFollow);
       ntrk++;
     }
     delete t;
@@ -1627,8 +1640,9 @@ AliCluster *AliITStrackerUpgrade::GetCluster(Int_t index) const {
   //       Return pointer to a given cluster
   //--------------------------------------------------------------------
   Int_t l=(index & 0xf0000000) >> 28;
+  Int_t c=(index & 0x0fffffff) >> 0;
   AliInfo(Form("index %i  cluster index %i layer %i", index,index & 0x0fffffff,index & 0xf0000000));
-  return fLayers[l]->GetCluster(index);
+  return fLayers[l]->GetCluster(c);
 }
 //______________________________________________________________________________
 Int_t AliITStrackerUpgrade::CorrectForPipeMaterial(AliITStrackMI *t, TString direction) {
@@ -1677,15 +1691,15 @@ if (!t->GetLocalXat(rToGo,xToGo)) return 0;
 Double_t xOverX0,x0,lengthTimesMeanDensity;
 
 switch(mode) {
- case 1:
+ case 0:
    xOverX0 = AliITSRecoParam::GetdPipe();
    x0 = AliITSRecoParam::GetX0Be();
    lengthTimesMeanDensity = xOverX0*x0;
    lengthTimesMeanDensity *= dir;
    if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
    break;
- case 0:
-   if (!t->PropagateToTGeo(xToGo,1)) return 0;
+ case 1:
+   if (!t->PropagateToTGeo(xToGo,20)) return 0; // "20" removes the d0 mean shift of ~20 um -> To be understood.   
    break;
  case 2:
    if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");