ITS/UPGRADE/AliITSUVertexer.cxx

   1 #include <Riostream.h>
   2 #include <TBranch.h>
   3 #include <TClonesArray.h>
   4 #include <TMath.h>
   5 #include <TString.h>
   6 #include <TTree.h>
   7 #include "AliESDVertex.h"
   8 #include "AliITSUClusterLines.h"
   9 #include "AliITSUClusterPix.h"
  10 #include "AliITSUVertexer.h"
  11 #include "AliLog.h"
  12 #include "AliStrLine.h"
  13 #include "AliVertexerTracks.h"
  14
  15 using TMath::Abs;
  16 using TMath::Sqrt;
  17 using TMath::ATan2;
  18 using TMath::TwoPi;
  19 using TMath::BubbleLow;
  20
  21 using std::cout;
  22 using std::endl;
  23
  24 //////////////////////////////////////////////////////////////////////
  25 // This class is used to compute the position of all the primary    //
  26 // vertices in a single event using the upgraded ITS.               //
  27 // Optimizations ongoing.                                           //
  28 // Origin puccio@to.infn.it  Feb. 20 2014                           //
  29 //////////////////////////////////////////////////////////////////////
  30
  31 //_____________________________________________________________________________________________
  32 AliITSUVertexer::AliITSUVertexer(Double_t phicut, Double_t zcut, Double_t paircut,Double_t clcut, Int_t cclcut) : AliVertexer(),
  33                                                                                                                   fClusterContribCut(cclcut),
  34                                                                                                                   fClusterCut(clcut),
  35                                                                                                                   fClusterIndex(),
  36                                                                                                                   fClusterPhi(),
  37                                                                                                                   fClusters(),
  38                                                                                                                   fLines("AliStrLine",1000),
  39                                                                                                                   fLinesClusters("AliITSUClusterLines.h",1000),
  40                                                                                                                   fLinesPhi(0),
  41                                                                                                                   fNoClusters(0),
  42                                                                                                                   fNoLines(0),
  43                                                                                                                   fNoVertices(0),
  44                                                                                                                   fPairCut(paircut),
  45                                                                                                                   fPhiCut(phicut),
  46                                                                                                                   fZCut(zcut),
  47                                                                                                                   fUsedClusters(),
  48                                                                                                                   fUsedLines(),
  49                                                                                                                   fVertices(NULL)
  50 #ifdef MC_CHECK
  51 ,fGoodLines(0),fGoodLinesPhi(0),fParticleId(0)
  52 #endif
  53 {
  54   // Standard I/O constructor
  55 }
  56
  57 //_____________________________________________________________________________________________
  58 AliITSUVertexer::~AliITSUVertexer() {
  59   // Destructor
  60   Reset();
  61 }
  62
  63 //_____________________________________________________________________________________________
  64 void AliITSUVertexer::FindVerticesForCurrentEvent() {
  65   // Try to find all the primary vertices in the current
  66   fNoVertices=0;
  67   FindTracklets();
  68   if(fNoLines<2) {
  69     //fVertices.push_back(AliESDVertex());
  70     return;// AliESDVertex();
  71   }
  72
  73   //  fVertices.push_back(AliVertexerTracks::TrackletVertexFinder(&fLines,1));
  74   //fNoVertices=1;
  75   fUsedLines=new Short_t[fNoLines];
  76   for(UInt_t i=0;i<fNoLines;++i) fUsedLines[i]=-1;
  77
  78   fNoClusters=0;
  79   for(UInt_t i1=0;i1<fNoLines;++i1) {
  80     if(fUsedLines[i1]!=-1) continue;
  81     AliStrLine* line1 = (AliStrLine*)fLines.At(i1);
  82     for(UInt_t i2=i1+1;i2<fNoLines;++i2) {
  83       if(fUsedLines[i2]!=-1) continue;
  84       AliStrLine* line2 = (AliStrLine*)fLines.At(i2);
  85       if(line1->GetDCA(line2)<=fPairCut) {
  86         //cout << fNoClusters <<" " << i1 << " " << i2 << " ";
  87         new(fLinesClusters[fNoClusters])AliITSUClusterLines(i1,line1,i2,line2);
  88         AliITSUClusterLines* current=(AliITSUClusterLines*)fLinesClusters.At(fNoClusters);
  89         Double_t p[3];
  90         current->GetVertex(p);
  91         if((p[0]*p[0]+p[1]*p[1])>=4) { // Beam pipe check
  92           fLinesClusters.RemoveAt(fNoClusters);
  93           fLinesClusters.Compress();
  94           break;
  95         }
  96         fUsedLines[i1]=fNoClusters;
  97         fUsedLines[i2]=fNoClusters;
  98         for(UInt_t i3=0;i3<fNoLines;++i3) {
  99           if(fUsedLines[i3]!=-1) continue;
 100           AliStrLine *line3 = (AliStrLine*)fLines.At(i3);
 101           //cout << p[0] << " " << p[1] << " " << p[2] << endl;
 102           //line3->PrintStatus();
 103           if(line3->GetDistFromPoint(p)<=fPairCut) {
 104             //cout << i3 << " ";
 105             current->Add(i3,line3);
 106             fUsedLines[i3]=fNoClusters;
 107             current->GetVertex(p);
 108           }
 109         }
 110         ++fNoClusters;
 111         //cout << endl;
 112         break;
 113       }
 114     }
 115   }
 116
 117   fLinesClusters.Sort();
 118
 119   for(UInt_t i0=0;i0<fNoClusters; ++i0) {
 120     Double_t p0[3],p1[3];
 121     AliITSUClusterLines *clu0 = (AliITSUClusterLines*)fLinesClusters.At(i0);
 122     clu0->GetVertex(p0);
 123     for(UInt_t i1=i0+1;i1<fNoClusters; ++i1) {
 124       AliITSUClusterLines *clu1 = (AliITSUClusterLines*)fLinesClusters.At(i1);
 125       clu1->GetVertex(p1);
 126       if (TMath::Abs(p0[2]-p1[2])<=fClusterCut) {
 127         Double_t distance=(p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1])+(p0[2]-p1[2])*(p0[2]-p1[2]);
 128         //Bool_t flag=kFALSE;
 129         if(distance<=fPairCut*fPairCut) {
 130           UInt_t n=0;
 131           Int_t *labels=clu1->GetLabels(n);
 132           for(UInt_t icl=0; icl<n; ++icl) clu0->Add(labels[icl],(AliStrLine*)fLines.At(labels[icl]));
 133           clu0->GetVertex(p0);
 134           //flag=kTRUE;
 135         }
 136         fLinesClusters.RemoveAt(i1);
 137         fLinesClusters.Compress();
 138         fNoClusters--;
 139         i1--;
 140         //if(flag) i1=10;
 141       }
 142     }
 143   }
 144
 145   fVertices=new AliESDVertex[fNoClusters];
 146   for(UInt_t i0=0; i0<fNoClusters; ++i0) {
 147     AliITSUClusterLines *clu0 = (AliITSUClusterLines*)fLinesClusters.At(i0);
 148     Int_t size=clu0->GetSize();
 149     if(size<fClusterContribCut&&fNoClusters>1) {
 150       fLinesClusters.RemoveAt(i0);
 151       fLinesClusters.Compress();
 152       fNoClusters--;
 153       continue;
 154     }
 155     Double_t p0[3],cov[6];
 156     clu0->GetVertex(p0);
 157     clu0->GetCovMatrix(cov);
 158     if((p0[0]*p0[0]+p0[1]*p0[1])<1.98*1.98) {
 159       fVertices[fNoVertices++]=AliESDVertex(p0,cov,99999.,size);
 160     }
 161   }
 162
 163   return;// AliVertexerTracks::TrackletVertexFinder(&fLines,0);
 164 }
 165
 166 //______________________________________________________________________
 167 AliESDVertex* AliITSUVertexer::FindVertexForCurrentEvent(TTree *cluTree)
 168 {
 169 // Reconstruction of all the primary vertices in the event. It returns the vertex with the highest number of contributors.
 170   Reset();
 171   for(Int_t i=0;i<3;++i) {
 172     TBranch* br = cluTree->GetBranch(Form("ITSRecPoints%d",i));
 173     if (!br) return NULL;
 174     br->SetAddress(&fClusters[i]);
 175   }
 176   cluTree->GetEntry(0);
 177
 178   SortClusters();
 179
 180   FindVerticesForCurrentEvent();
 181   if(fNoVertices<1) return NULL;
 182   return new AliESDVertex(fVertices[0]);
 183 }
 184
 185 //_____________________________________________________________________________________________
 186 Int_t AliITSUVertexer::MatchPoints(UShort_t layer, Double_t anchor, Double_t *p0, Double_t *p1) {
 187   // Method for matching clusters with similar phi and z inside the range [zmin,zmax]
 188   AliDebug(3,Form("Matching points on layer %i...\n",layer));
 189
 190   Double_t xyz[3][3];        // {{0: 'min', 1: 'mid', 2: 'max'},{0: 'x', 1: 'y', 2: 'z'}}
 191   AliITSUClusterPix* cl[3];  // {0: 'min', 1: 'mid', 2: 'max'}
 192   Double_t phi[3];           // {0: 'min', 1: 'mid', 2: 'max'}
 193
 194   Int_t nocl=fClusters[layer]->GetEntriesFast();
 195   Int_t imin=0,imax=nocl!=0 ? nocl-1 : nocl,imid=(imax+imin)/2;
 196
 197   Double_t a=0,b=0,c=0;
 198   if(layer==2) {
 199     a=(p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1]);
 200     b=(p1[0]-p0[0])*p0[0]+(p1[1]-p0[1])*p0[1];
 201     c=p0[0]*p0[0]+p0[1]*p0[1];
 202   }
 203
 204   Int_t flag=-1;
 205
 206   while (imax > imin) {
 207     while(fUsedClusters[layer][fClusterIndex[layer][imin]]==kTRUE&&imin<imax) ++imin;
 208     while(fUsedClusters[layer][fClusterIndex[layer][imax]]==kTRUE&&imax>0) --imax;
 209     if(imax<imin) return -1;
 210
 211     cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imin]);
 212     cl[0]->GetGlobalXYZ(xyz[0]);
 213     phi[0] = fClusterPhi[layer][fClusterIndex[layer][imin]];
 214     AliDebug(4,Form("Cluster min: %i, %i, %f, %i\n",imin,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][imin]]));
 215     if((Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) {
 216       AliDebug(4,Form("Ok min: %i \n",imin));
 217       if(layer==2) {
 218         if(fUsedClusters[layer][fClusterIndex[layer][imin]]) {
 219           flag=1;
 220           break;
 221         }
 222         c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
 223         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
 224         Double_t zmin=z-fZCut,zmax=z+fZCut;
 225         if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
 226           AliDebug(4,Form("Ok Z: %i \n",imin));
 227           return fClusterIndex[layer][imin];
 228         }
 229         AliDebug(4,Form("No Z match: %i \n",imin));
 230         c=p0[0]*p0[0]+p0[1]*p0[1];
 231         flag=1;
 232         break;
 233       } else if(!fUsedClusters[layer][fClusterIndex[layer][imin]]) return fClusterIndex[layer][imin];
 234     }
 235
 236     cl[2] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imax]);
 237     cl[2]->GetGlobalXYZ(xyz[2]);
 238     phi[2] = fClusterPhi[layer][fClusterIndex[layer][imax]];
 239     AliDebug(4,Form("Cluster max: %i, %i, %f, %i\n",imax,cl[2]->GetLabel(0),phi[2]-anchor,fUsedClusters[layer][fClusterIndex[layer][imax]]));
 240     if((Abs(phi[2]-anchor)<=fPhiCut||Abs(Abs(phi[2]-anchor)-TMath::TwoPi())<=fPhiCut)) {
 241       AliDebug(4,Form("Ok max: %i \n",imax));
 242       if(layer==2) {
 243         if(fUsedClusters[layer][fClusterIndex[layer][imax]]) {
 244           flag=3;
 245           break;
 246         }
 247         c-=(xyz[2][0]*xyz[2][0]+xyz[2][1]*xyz[2][1]);
 248         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
 249         Double_t zmin=z-fZCut,zmax=z+fZCut;
 250         if(zmin<=xyz[2][2]&&xyz[2][2]<=zmax) {
 251           AliDebug(4,Form("Ok Z: %i \n",imax));
 252           return fClusterIndex[layer][imax];
 253         }
 254         c=p0[0]*p0[0]+p0[1]*p0[1];
 255         AliDebug(4,Form("No Z match: %i \n",imax));
 256         flag=3;
 257         break;
 258       } else if(!fUsedClusters[layer][fClusterIndex[layer][imax]]) return fClusterIndex[layer][imax];
 259     }
 260
 261     imid=(imax+imin)/2;
 262     if(imid==imin) return -1;
 263     Int_t step=1,sign=-1;
 264     while(fUsedClusters[layer][fClusterIndex[layer][imid]]) {
 265       imid=imid+step;
 266       sign*=-1;
 267       step=(step+sign)*(-1);
 268       if(imid>=imax||imid<=imin) return -1;
 269     }
 270
 271     cl[1] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imid]);
 272     cl[1]->GetGlobalXYZ(xyz[1]);
 273     phi[1] = fClusterPhi[layer][fClusterIndex[layer][imid]];
 274     AliDebug(4,Form("Cluster mid: %i, %i, %f, %i \n",imid,cl[1]->GetLabel(0),phi[1]-anchor,fUsedClusters[layer][fClusterIndex[layer][imid]]));
 275     //cout << imin << " " << imid << " " << imax << endl;
 276     //cout << fClusterIndex[layer][imid] << endl;
 277     if((Abs(phi[1]-anchor)<=fPhiCut||Abs(Abs(phi[1]-anchor)-TMath::TwoPi())<=fPhiCut)) {
 278       AliDebug(4,Form("Ok mid: %i \n",imid));
 279       if(layer==2) {
 280         c-=(xyz[1][0]*xyz[1][0]+xyz[1][1]*xyz[1][1]);
 281         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
 282         Double_t zmin=z-fZCut,zmax=z+fZCut;
 283         if(zmin<=xyz[1][2]&&xyz[1][2]<=zmax) {
 284           AliDebug(4,Form("Ok Z: %i \n",imid));
 285           return fClusterIndex[layer][imid];
 286         }
 287         AliDebug(4,Form("No Z match: %i \n",imid));
 288         c=p0[0]*p0[0]+p0[1]*p0[1];
 289         flag=2;
 290         break;
 291       } else if(!fUsedClusters[layer][fClusterIndex[layer][imid]]) return fClusterIndex[layer][imid];
 292     }
 293
 294     // determine which subarray to search
 295     if (phi[1] < anchor) {
 296       // change min index to search upper subarray
 297       AliDebug(4,"Case minor\n");
 298       imin = imid + 1;
 299       --imax;
 300     } else if (phi[1] > anchor) {
 301       AliDebug(4,"Case greater\n");
 302       // change max index to search lower subarray
 303       ++imin;
 304       imax = imid - 1;
 305     } else if(!fUsedClusters[layer][fClusterIndex[layer][imid]]) {
 306       return fClusterIndex[layer][imid];
 307     } else return -1;
 308   }
 309
 310   if(flag>-1) {
 311     AliDebug(4,"Flag issued, starting forward backward check\n");
 312     Int_t start=imid;
 313     switch(flag) {
 314     case 1:
 315       start=imin;
 316       break;
 317     case 2:
 318       start=imid;
 319       break;
 320     case 3:
 321       start=imax;
 322       break;
 323     }
 324
 325     Int_t curr=start-1;
 326     Bool_t lap=kFALSE;
 327     while(1){
 328       if(curr==-1&&!lap) {
 329         lap=kTRUE;
 330         curr=nocl-1;
 331       } else if(lap) break;
 332       cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][curr]);
 333       cl[0]->GetGlobalXYZ(xyz[0]);
 334       phi[0] = fClusterPhi[layer][fClusterIndex[layer][curr]];
 335       AliDebug(4,Form("Looking backward: %i, %i, %f, %i\n",curr,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][curr]]));
 336       if(!(Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) break;
 337       if(fUsedClusters[layer][fClusterIndex[layer][curr]]) {
 338         curr--;
 339         continue;
 340       }
 341
 342       AliDebug(4,Form("Ok backward: %i \n",curr));
 343
 344       c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
 345       Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
 346       Double_t zmin=z-fZCut,zmax=z+fZCut;
 347       if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
 348         AliDebug(4,Form("Ok Z: %i \n",curr));
 349         return fClusterIndex[layer][curr];
 350       }
 351       AliDebug(4,Form("No Z match: %i \n",curr));
 352       c=p0[0]*p0[0]+p0[1]*p0[1];
 353       curr--;
 354     }
 355
 356     lap=kFALSE;
 357     curr=start+1;
 358     while(1){
 359       if(curr==nocl&&!lap) {
 360         curr=0;
 361         lap=kTRUE;
 362       } else if(lap) break;
 363       cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][curr]);
 364       cl[0]->GetGlobalXYZ(xyz[0]);
 365       phi[0] = fClusterPhi[layer][fClusterIndex[layer][curr]];
 366       AliDebug(4,Form("Looking forward: %i, %i, %f, %i\n",curr,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][curr]]));
 367       if(!(Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)) break;
 368       if(fUsedClusters[layer][fClusterIndex[layer][curr]]) {
 369         curr++;
 370         continue;
 371       }
 372       AliDebug(4,Form("Ok forward: %i \n",curr));
 373       c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
 374       Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
 375       Double_t zmin=z-fZCut,zmax=z+fZCut;
 376       if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
 377         AliDebug(4,Form("Ok Z: %i \n",curr));
 378         return fClusterIndex[layer][curr];
 379       }
 380       AliDebug(4,Form("No Z match: %i \n",curr));
 381       c=p0[0]*p0[0]+p0[1]*p0[1];
 382       curr++;
 383     }
 384   }
 385
 386   if(imax==imin&&imax!=0) {
 387     cl[0] = (AliITSUClusterPix*)fClusters[layer]->At(fClusterIndex[layer][imin]);
 388     cl[0]->GetGlobalXYZ(xyz[0]);
 389     phi[0] = fClusterPhi[layer][fClusterIndex[layer][imin]];
 390     AliDebug(4,Form("Cluster eq: %i, %i, %f, %i\n",imin,cl[0]->GetLabel(0),phi[0]-anchor,fUsedClusters[layer][fClusterIndex[layer][imin]]));
 391     if((Abs(phi[0]-anchor)<=fPhiCut||Abs(Abs(phi[0]-anchor)-TMath::TwoPi())<=fPhiCut)&&!fUsedClusters[layer][fClusterIndex[layer][imin]]) {
 392       AliDebug(4,Form("Ok eq: %i \n",imin));
 393       if(layer==2) {
 394         c-=(xyz[0][0]*xyz[0][0]+xyz[0][1]*xyz[0][1]);
 395         Double_t z=p0[2]+(p1[2]-p0[2])*(-b+Sqrt(b*b-a*c))/a;
 396         Double_t zmin=z-fZCut,zmax=z+fZCut;
 397         if(zmin<=xyz[0][2]&&xyz[0][2]<=zmax) {
 398           AliDebug(4,Form("Ok Z eq: %i \n",imin));
 399           return fClusterIndex[layer][imin];
 400         }
 401         AliDebug(4,Form("No Z eq: %i \n",imin));
 402         c=p0[0]*p0[0]+p0[1]*p0[1];
 403       } else return fClusterIndex[layer][imin];
 404     }
 405   }
 406   // no match found :(
 407   return -1;
 408
 409 }
 410
 411 //_____________________________________________________________________________________________
 412 void AliITSUVertexer::PrintStatus() const {
 413   // Prints all the cuts and important data members for the current status
 414   cout << "Cut on phi: " << fPhiCut << endl;
 415   cout << "Cut on z: " << fZCut << endl;
 416 }
 417
 418 //_____________________________________________________________________________________________
 419 void AliITSUVertexer::Reset() {
 420   // Resets the vertexer for a new event (or for its destruction)
 421   AliDebug(2,"Resetting the vertexer...\n");
 422   fNoVertices=0;
 423   for(Int_t i=0;i<3;++i) {
 424     delete[] fUsedClusters[i];
 425     delete[] fClusterIndex[i];
 426     delete[] fClusterPhi[i];
 427   }
 428   if(fNoLines>2) {
 429     delete []fUsedLines;
 430   }
 431
 432   delete[] fVertices;
 433
 434   fLinesPhi=0;
 435   fLines.Clear();
 436   fLinesClusters.Clear();
 437
 438   #ifdef MC_CHECK
 439   fGoodLines=0;
 440   fGoodLinesPhi=0;
 441   delete[] fParticleId;
 442   #endif
 443
 444   //delete fVertices;
 445 }
 446
 447 //_____________________________________________________________________________________________
 448 void AliITSUVertexer::FindTracklets() {
 449   // It combines recpoints over the first three layers to define a list of tracklets
 450   // Here it uses the ordered list of points. Second and third layers are ordered using phi
 451   AliDebug(2,"Calling the trackleter...\n");
 452
 453   UInt_t noPntL[3];
 454   for(UShort_t i=0;i<3;++i) {
 455     noPntL[i]=fClusters[i]->GetEntries();
 456     fUsedClusters[i]=new Bool_t[noPntL[i]];
 457     for(UInt_t ii=0;ii<noPntL[i];++ii) fUsedClusters[i][ii]=kFALSE;
 458   }
 459
 460   #ifdef MC_CHECK
 461   fParticleId=new UInt_t[noPntL[0]];
 462   #endif
 463
 464   fNoLines=0;
 465   UInt_t nolinesphi=0;
 466   Double_t p0[3],p1[3],pp2[3];
 467   for(UInt_t i0=0;i0<noPntL[0];++i0) {
 468     if(fUsedClusters[0][i0]) continue;
 469     AliITSUClusterPix* cluster0=(AliITSUClusterPix*)fClusters[0]->At(i0);
 470     cluster0->GetGlobalXYZ(p0);
 471     vector<Int_t> tmp;
 472     Int_t i1=0;
 473     Int_t label0=cluster0->GetLabel(0);
 474     AliDebug(4,Form("Layer 0: %i\n",label0));
 475     while(i1>=0) {
 476       i1 = MatchPoints(1,fClusterPhi[0][i0]);
 477       if(i1<0) break;
 478
 479       ++nolinesphi;
 480       AliITSUClusterPix* cluster1=(AliITSUClusterPix*)fClusters[1]->At(i1);
 481       cluster1->GetGlobalXYZ(p1);
 482       tmp.push_back(i1);
 483       fUsedClusters[1][i1]=kTRUE;
 484
 485       Int_t i2 = MatchPoints(2,fClusterPhi[1][i1],p0,p1);
 486       if(i2<0) continue;
 487
 488       #ifdef MC_CHECK
 489       CheckMC(i0,i1,i2);
 490       #endif
 491
 492       AliITSUClusterPix* cluster2=(AliITSUClusterPix*)fClusters[2]->At(i2);
 493       cluster2->GetGlobalXYZ(pp2);
 494
 495       fUsedClusters[0][i0]=kTRUE;
 496       fUsedClusters[2][i2]=kTRUE;
 497
 498       // Errors to be checked...
 499
 500       Float_t cov0[6],cov1[6],cov2[6];
 501       cluster0->GetGlobalCov(cov0);
 502       cluster1->GetGlobalCov(cov1);
 503       cluster2->GetGlobalCov(cov2);
 504
 505       //Error on tracklet direction near the vertex
 506       Double_t rad1=TMath::Sqrt(p0[0]*p0[0]+p0[1]*p0[1]);
 507       Double_t rad2=TMath::Sqrt(p1[0]*p1[0]+p1[1]*p1[1]);
 508       Double_t factor=(rad1+rad2)/(rad2-rad1);
 509
 510       //Curvature error
 511       Double_t curvErr=0;
 512       Double_t bField=0.5;
 513       Double_t meanPtSelTrk=0.630;
 514       Double_t curvRadius=meanPtSelTrk/(0.3*bField)*100; //cm
 515       Double_t dRad=TMath::Sqrt((p0[0]-p1[0])*(p0[0]-p1[0])+(p0[1]-p1[1])*(p0[1]-p1[1]));
 516       Double_t aux=dRad/2.+rad1;
 517       curvErr=TMath::Sqrt(curvRadius*curvRadius-dRad*dRad/4.)-TMath::Sqrt(curvRadius*curvRadius-aux*aux); //cm
 518
 519       Double_t sq[3],wmat[9]={1,0,0,0,1,0,0,0,1};
 520       sq[0]=(cov0[0]+curvErr*curvErr/2.)*factor*factor;//cov1[0]+cov2[0]);
 521       sq[1]=(cov0[3]+curvErr*curvErr/2.)*factor*factor;//cov1[3]+cov2[3]);
 522       sq[2]=(cov0[5])*factor*factor;//cov1[5]+cov2[5]);
 523
 524       // Multiple scattering
 525       Double_t meanPSelTrk=0.875;
 526       Double_t pOverMass=meanPSelTrk/0.140;
 527       Double_t beta2=pOverMass*pOverMass/(1+pOverMass*pOverMass);
 528       Double_t p2=meanPSelTrk*meanPSelTrk;
 529       Double_t rBP=1.98; // Beam pipe radius
 530       Double_t dBP=0.08/35.3; // 800 um of Be
 531       Double_t dL1=0.01; //approx. 1% of radiation length
 532       Double_t theta2BP=14.1*14.1/(beta2*p2*1e6)*dBP;
 533       Double_t theta2L1=14.1*14.1/(beta2*p2*1e6)*dL1;
 534       Double_t rtantheta1=(rad2-rad1)*TMath::Tan(TMath::Sqrt(theta2L1));
 535       Double_t rtanthetaBP=(rad1-rBP)*TMath::Tan(TMath::Sqrt(theta2BP));
 536       for(Int_t ico=0; ico<3;ico++){
 537         sq[ico]+=rtantheta1*rtantheta1*factor*factor/3.;
 538         sq[ico]+=rtanthetaBP*rtanthetaBP*factor*factor/3.;
 539       }
 540
 541       if(sq[0]!=0) wmat[0]=1/sq[0];
 542       if(sq[1]!=0) wmat[4]=1/sq[1];
 543       if(sq[2]!=0) wmat[8]=1/sq[2];
 544
 545       /*Int_t label0=cluster0->GetLabel(0);
 546       Int_t label1=cluster1->GetLabel(0);
 547       Int_t label2=cluster2->GetLabel(0);*/
 548       //cout << label0 << " " << label1 << " "<<label2 << endl;
 549
 550       //printf("%f\t%f\t%f\n-\t%f\t%f\n-\t-\t%f\n",wmat[0],wmat[1],wmat[2],wmat[3],wmat[4],wmat[5]);
 551       //gSystem->Exec(Form("echo %i >> trkl",cluster0->GetLabel(0)+cluster1->GetLabel(0)));
 552
 553       new(fLines[fNoLines++])AliStrLine(p0,sq,wmat,p1,kTRUE);
 554       break;
 555     }
 556     for(UInt_t itmp=0; itmp<tmp.size(); ++itmp) fUsedClusters[1][tmp.at(itmp)]=kFALSE;
 557     tmp.clear();
 558     //((AliStrLine*)fLines[fNoLines-1])->PrintStatus();
 559     //printf("(%f,%f,%f) and (%f,%f,%f)\n",p0[0],p0[1],p0[2],p1[0],p1[1],p1[2]);
 560   }
 561
 562   fLinesPhi=nolinesphi;
 563   //cout << Form("=======================================================\nNolinesphi: %i, nolines: %i\nGood Nolinesphi: %i, good nolines: %i\n",nolinesphi,fNoLines,fGoodLinesPhi,fGoodLines);
 564 }
 565
 566 //___________________________________________________________________________
 567 void AliITSUVertexer::SortClusters() {
 568   // Reading of the clusters on the first three layer of the upgraded ITS
 569   for(Int_t i=0;i<3;++i) {
 570     TClonesArray *clr=fClusters[i];
 571     Int_t nocl=clr->GetEntriesFast();
 572     if(nocl==0) {
 573       fClusterPhi[i]=new Double_t[1];
 574       fClusterPhi[i][0]=-999999;
 575       fClusterIndex[i]=new Int_t[1];
 576       fClusterIndex[i][0]=0;
 577     } else {
 578       fClusterPhi[i]=new Double_t[nocl];
 579       fClusterIndex[i]=new Int_t[nocl];
 580       for(Int_t k=0;k<nocl;++k) {
 581         AliITSUClusterPix* cl=(AliITSUClusterPix*)clr->At(k);
 582         Double_t pt[3];
 583         cl->GetGlobalXYZ(pt);
 584         fClusterPhi[i][k]=ATan2(pt[1],pt[0]);
 585       }
 586       BubbleLow(nocl,fClusterPhi[i],fClusterIndex[i]);
 587     }
 588   }
 589 }
 590
 591 #ifdef MC_CHECK
 592 //___________________________________________________________________________
 593 Bool_t AliITSUVertexer::CheckMC(UInt_t i0, UInt_t i1, UInt_t i2) {
 594   // Debugging function
 595   //cout << "Checking MC truth" << endl;
 596   Int_t label=0;
 597   Bool_t flag=kFALSE;
 598   AliITSUClusterPix* p0=(AliITSUClusterPix*)fClusters[0]->At(i0);
 599   AliITSUClusterPix* p1=(AliITSUClusterPix*)fClusters[1]->At(i1);
 600   for(Int_t i=0; i<3; ++i) {
 601     label=p0->GetLabel(i);
 602     for(Int_t j=0; j<3; ++j) {
 603       if(label==p1->GetLabel(j)&&label>0) {
 604         fGoodLinesPhi++;
 605         flag=kTRUE;
 606         break;
 607       }
 608     }
 609     if(flag) break;
 610   }
 611
 612   if(!flag) return kFALSE;
 613
 614   AliITSUClusterPix* p2=(AliITSUClusterPix*)fClusters[2]->At(i2);
 615   for(Int_t j=0; j<3; ++j) {
 616     if(label==p2->GetLabel(j)) {
 617       fParticleId[fGoodLines]=label;
 618       fGoodLines++;
 619       //cout << label << endl;
 620       return kTRUE;
 621     }
 622   }
 623   return kFALSE;
 624 }
 625 #endif