AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliMFTTrackFinder.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 //============================================================================================
17 //
18 // Class for finding Tracks from Cluster of the ALICE Muon Forward Tracker
19 //
20 // Contact author: raphael.tieulent@cern.ch
21 //
22 //============================================================================================
23 
24 #include <fstream>
25 #include <iostream>
26 
27 #include "TMath.h"
28 #include "TF1.h"
29 
30 #include "AliCodeTimer.h"
31 #include "AliLog.h"
32 #include "AliMFTCAHit.h"
33 #include "AliMFTCARoad.h"
34 #include "AliMFTCACell.h"
35 #include "AliMFTCALayer.h"
36 #include "AliMFTCATrack.h"
37 #include "AliMFTCluster.h"
38 #include "AliMFTTrackFinder.h"
39 
40 
42 
43 //___________________________________________________________________________
45 TObject(),
46 fXCut(0.0025),
47 fYCut(0.0025),
48 fMaxSegAngle(180.-15.),
49 fCellGID(0),
50 fMaxCellStatus(0),
51 fNlayers(fNDetMax),
52 fNtracks(0),
53 fCalcVertex(kFALSE),
54 fZVertCalc(0.),
55 fZVertDet(0.),
56 fZVertRange(),
57 fMinTrackLength(fNDetMax),
58 fPlaneDetEff(),
59 fAddNoise(kTRUE),
60 fPixelNoise(1.E-5),
61 fAddQED(kFALSE),
62 fMBRate(0.),
63 fCMOSIntTime(0.),
64 fThick(0.004),
65 fReadGeom(kTRUE),
66 fUseTF(kFALSE),
67 fDebug(0),
68 fCPUTime(0.),
69 fRealTime(0.),
70 fNDifTracks(0),
71 fPlanesZ(),
72 fZGap(1.0),
73 fNRoads(0),
74 fErrX(0.),
75 fErrY(0.)
76 {
77 
78 #ifdef OLDGEOM
79  fZGap = 0.2;
80 #else
81  fZGap = 1.0;
82 #endif
83 
84  fZVertRange[0] = fZVertRange[1] = 0.;
85 
86  fHistList = new TList();
87 
88  fTracks = new TClonesArray("AliMFTCATrack", 1000);
89 
90  fRoads = new TClonesArray("AliMFTCARoad", 1000);
91 
92  for (Int_t iL = 0; iL < fNlayers; iL++) {
93 
94  fACutV[iL] = 0.20;
95  fACutN[iL] = 0.10;
96 
97  fLayers[iL] = new AliMFTCALayer();
98  fLayers[iL]->SetID(iL);
99  fPlaneDetEff[iL] = 1.00;
100  fPlanesZ[iL] = 0.;
101 
102  hDA[iL] = new TH1F(Form("hDA%d",iL),Form("#Delta #theta between two segments in Layer %d; #Delta #theta (deg)",iL),200,0.,0.5);
103  hDAv[iL] = new TH1F(Form("hDAv%d",iL),Form("#Delta #theta with respect to the vertex in Layer %d; #Delta #theta (deg)",iL),200,0.,3.0);
104  fHistList->Add(hDA[iL]);
105  fHistList->Add(hDAv[iL]);
106  hDXY[iL] = new TH2F(Form("hDXY%d",iL),Form("Distance between 2 hits on layer %d ; X (microns) ; Y (microns)",iL),100,0.,10.,100,0.,+10.);
107  //hDXY[iL] = new TH2F(Form("hDXY%d",iL),Form("Distance between 2 hits on layer %d ; X (cm) ; Y (cm)",iL),100,0.,+1.0,100,0.,+1.0);
108  fHistList->Add(hDXY[iL]);
109 
110  }
111  hNGoodCell = new TH1F("hNGoodCell","Number of good cell in the track",5,-0.5,4.5);
112  fHistList->Add(hNGoodCell);
113  hTrackType = new TH1F("hTrackType","Track type",4,-0.5,3.5);
114  fHistList->Add(hTrackType);
115  hAngleCells = new TH1F("hAngleCell1gap2gaps","#Delta #theta between two segments one 1 gap and one 2 layers gap; #Delta #theta (deg)",200,0.,45.);
116  fHistList->Add(hAngleCells);
117  hThetaCells = new TH1F("hThetaCells","#theta of the cells; #theta (deg)",200,0.,15.);
118  fHistList->Add(hThetaCells);
119 
120  // limited range for possible z vertex
121  hVertZ = new TH1F("hVertZ","hVertZ",400,-10.,+10.);
122  hVertZa = new TH1F("hVertZa","hVertZa",400,-10.,+10.);
123 
124  fHistList->Add(hVertZ);
125  fHistList->Add(hVertZa);
126 
127  // temporary
128  hHitDifXY = new TH2F("hHitDifXY","hHitDifXY",100,-0.04,+0.04,100,-0.04,+0.04);
129  fHistList->Add(hHitDifXY);
130  hAngDifAll = new TH1F("hAngDifAll","hAngDifAll",1000,0.,+0.5);
131  fHistList->Add(hAngDifAll);
132  hAngDifDup = new TH1F("hAngDifDup","hAngDifDup",1000,0.,+0.5);
133  fHistList->Add(hAngDifDup);
134  hIntDifXYAll = new TH2F("hIntDifXYAll","hIntDifXYAll",100,-0.01,+0.01,100,-0.01,+0.01);
135  fHistList->Add(hIntDifXYAll);
136 #ifdef OLDGEOM
137  hIntDifXYDup = new TH2F("hIntDifXYDup","hIntDifXYDup",100,-0.01,+0.01,100,-0.01,+0.01);
138 #else
139  hIntDifXYDup = new TH2F("hIntDifXYDup","hIntDifXYDup",100,-0.1,+0.1,100,-0.1,+0.1);
140 #endif
141  fHistList->Add(hIntDifXYDup);
142 
143 }
144 
145 //___________________________________________________________________________
147 {
148 
149 }
150 //___________________________________________________________________________
151 void AliMFTTrackFinder::Init(Char_t *parfile)
152 {
153 
154  AliInfo(Form("Parameter file set to %s ",parfile));
155 
156  ReadParam(parfile);
157 
158  ReadGeom();
159 
160  SetDebug(1);
161 
162 
163 
164 }
165 //___________________________________________________________________________
166 void AliMFTTrackFinder::LoadClusters( TClonesArray *clusterArrayFront[AliMFTConstants::fNMaxPlanes], TClonesArray *clusterArrayBack[AliMFTConstants::fNMaxPlanes] ){
167  AliCodeTimerAuto("",0);
168  AliMFTCALayer *caLayer = NULL;
169  AliMFTCAHit *caHit = NULL;
170 
171  for (int iPlane = 0; iPlane<AliMFTConstants::kNDisks; iPlane++) {
172  Int_t nClusterFront = clusterArrayFront[iPlane]->GetEntriesFast();
173  Int_t nClusterBack = clusterArrayBack[iPlane]->GetEntriesFast();
174  AliDebug(1,Form("Loading %d + %d = %d Clusters for Plane %d",nClusterFront , nClusterBack,nClusterFront + nClusterBack, iPlane));
175 
176  // Treating FRONT face of the plane
177  caLayer = GetLayer(iPlane*2);
178  for (Int_t iCluster=0; iCluster<nClusterFront; iCluster++) {
179  caHit = caLayer->AddHit();
180  AliMFTCluster * cluster = (AliMFTCluster *)clusterArrayFront[iPlane]->At(iCluster);
181  caHit->SetPos(cluster->GetX(), cluster->GetY(), cluster->GetZ());
182  caHit->SetTrackGID(cluster->GetMCLabel(0),iPlane*2,caLayer->GetNhits()-1,0);
183  caHit->SetMFTClsId(iCluster);
184  }
185  // Treating BACK face of the plane
186  caLayer = GetLayer(iPlane*2+1);
187  for (Int_t iCluster=0; iCluster<nClusterBack; iCluster++) {
188  caHit = caLayer->AddHit();
189  AliMFTCluster * cluster = (AliMFTCluster *)clusterArrayBack[iPlane]->At(iCluster);
190  caHit->SetPos(cluster->GetX(), cluster->GetY(), cluster->GetZ());
191  caHit->SetTrackGID(cluster->GetMCLabel(0),iPlane*2+1,caLayer->GetNhits()-1,0);
192  caHit->SetMFTClsId(iCluster);
193  }
194 
195  }
196 }
197 
198 //___________________________________________________________________________
199 void AliMFTTrackFinder::ReadParam(Char_t *parfile)
200 {
201 
202  AliInfo(Form("Reading Parameter File %s ",parfile));
203 
204  std::ifstream in;
205  in.open(parfile,std::ios::in);
206 
207  std::string line;
208 
209  getline(in,line);
210  in >> fUseTF;
211  printf("Use the TrackFinder: %d \n",fUseTF);
212 
213  getline(in,line); getline(in,line);
214  in >> fNlayers;
215  printf("Number of detecting planes: %d \n",fNlayers);
216 
217  getline(in,line); getline(in,line);
218  in >> fThick;
219  printf("Layer thickness in X0: %5.3f \n",fThick);
220 
221  getline(in,line); getline(in,line);
222  in >> fPixelNoise;
223  printf("Pixel noise: %4.2e \n",fPixelNoise);
224  in >> fAddNoise;
225  printf("Add noise: %d \n",fAddNoise);
226 
227  getline(in,line); getline(in,line);
228  in >> fMinTrackLength;
229  printf("fMinTrackLength: %d \n",fMinTrackLength);
230 
231  getline(in,line); getline(in,line);
232  in >> fXCut;
233  printf("fXCut: %6.4f cm\n",fXCut);
234  in >> fYCut;
235  printf("fYCut: %6.4f cm\n",fYCut);
236 
237  getline(in,line); getline(in,line);
238  in >> fMaxSegAngle;
239  printf("fMaxSegAngle: %4.1f deg\n",fMaxSegAngle);
240  fMaxSegAngle = 180. - fMaxSegAngle;
241 
242  getline(in,line); getline(in,line);
243  for (Int_t i = 0; i < fNlayers; i++) {
244  in >> fACutV[i];
245  printf("fACutV[%d]: %4.2f \n",i,fACutV[i]);
246  }
247 
248  getline(in,line); getline(in,line);
249  for (Int_t i = 0; i < fNlayers; i++) {
250  in >> fACutN[i];
251  printf("fACutN[%d]: %4.2f \n",i,fACutN[i]);
252  }
253 
254  getline(in,line); getline(in,line);
255  for (Int_t i = 0; i < fNlayers; i++) {
256  in >> fPlaneDetEff[i];
257  printf("PlaneDetEff[%d]: %4.2f \n",i,fPlaneDetEff[i]);
258  }
259 
260  getline(in,line); getline(in,line);
261  in >> fCalcVertex;
262  printf("Calculate vertex: %d \n",fCalcVertex);
263 
264  getline(in,line); getline(in,line);
265  in >> fAddQED;
266  printf("Add QED hits: %d \n",fAddQED);
267  in >> fMBRate;
268  printf("Hadronic MB Rate: %5.1f kHz\n",fMBRate);
269  in >> fCMOSIntTime;
270  printf("CMOS integration time: %5.1f microsec \n",fCMOSIntTime);
271 
272  getline(in,line); getline(in,line);
273  in >> fReadGeom;
274  printf("Read geometry: %d \n",fReadGeom);
275  in >> fGeomName;
276  if (fReadGeom) printf("... from file: %s \n",fGeomName.Data());
277 
278  in.close();
279  printf("==================================\n");
280 
281 }
282 
283 //___________________________________________________________________________
284 void AliMFTTrackFinder::Clear(Option_t *)
285 {
286 
287  if (fTracks) fTracks->Clear("C");
288 
289  for (Int_t iL = 0; iL < fNlayers; iL++) {
290  fLayers[iL]->Clear("");
291  }
292 
293  fCellGID = 0;
294  fMaxCellStatus = 0;
295  fNtracks = 0;
296 
297  if (fRoads) fRoads->Clear("C");
298 
299  fNRoads = 0;
300 
301 }
302 
303 //___________________________________________________________________________
305 {
306 
307  for (Int_t iL = 0; iL < fNlayers; iL++) {
308  fLayers[iL]->ClearCells();
309  }
310 
311  fCellGID = 0;
312  fMaxCellStatus = 0;
313 
314 }
315 
316 //___________________________________________________________________________
318 {
319 
320  AliMFTCALayer *layer;
321  AliMFTCACell *cell;
322 
323  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
324  layer = GetLayer(iL);
325  for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
326  cell = layer->GetCell(iC);
327  cell->Reset();
328  }
329  }
330 
331 }
332 
333 //___________________________________________________________________________
335 
336  Double_t *c1h1 = cell1->GetHit1();
337  Double_t *c1h2 = cell1->GetHit2();
338  Double_t *c2h1 = cell2->GetHit1();
339  Double_t *c2h2 = cell2->GetHit2();
340 
341  TVector3 seg1 = TVector3(c1h2[0]-c1h1[0],c1h2[1]-c1h1[1],c1h2[2]-c1h1[2]);
342  TVector3 seg2 = TVector3(c2h2[0]-c2h1[0],c2h2[1]-c2h1[1],c2h2[2]-c2h1[2]);
343 
344  return (seg1.Angle(seg2))*TMath::RadToDeg();
345 
346 }
347 
348 //___________________________________________________________________________
350 
351  Double_t *ch1 = cell->GetHit1();
352  Double_t *ch2 = cell->GetHit2();
353 
354  return ch1[0]+(ch2[0]-ch1[0])/(ch2[2]-ch1[2])*(z-ch1[2]);
355 
356 }
357 
358 //___________________________________________________________________________
360 
361  Double_t *ch1 = cell->GetHit1();
362  Double_t *ch2 = cell->GetHit2();
363 
364  return ch1[1]+(ch2[1]-ch1[1])/(ch2[2]-ch1[2])*(z-ch1[2]);
365 
366 }
367 
368 //___________________________________________________________________________
370 
371  printf("Total number of cells = %d \n",GetNcells());
372 
373  // analyze cell for aliroot input
374 
375  Double_t cellAngDifCut; // [deg]
376  Double_t cellIntDifCut; // [cm]
377 
378  cellAngDifCut = 0.03;
379  cellIntDifCut = 0.01;
380  //cellIntDifCut = fXCut; // old
381 
382  AliMFTCALayer *layer;
383  AliMFTCACell *cell1, *cell2, *cellM, *cell;
384  Int_t trkid1h1, trkid1h2, trkid2h1, trkid2h2;
385  Double_t xc1h1, xc1h2, xc2h1, xc2h2;
386  Double_t yc1h1, yc1h2, yc2h1, yc2h2;
387  Double_t zc1h1, zc1h2, zc2h1, zc2h2, zMin, zMax;
388  Double_t cellAngDif, cellIntDifX, cellIntDifY, zint;
389  Double_t cellDifX1, cellDifX2, cellDifY1, cellDifY2;
390  Double_t hit1[3], hit2[3];
391  Bool_t multCell1;
392  TList *multCell = new TList();
393 
394  const Int_t nMaxh = 100;
395  Double_t xTr[nMaxh], yTr[nMaxh], zTr[nMaxh];
396  Double_t ax, axe, bx, bxe, ay, aye, by, bye;
397  Double_t xTrErrDet = 0.0025/TMath::Sqrt(12.);
398  Double_t yTrErrDet = 0.0025/TMath::Sqrt(12.);
399  Double_t xTrErrMS = 0.00055; // estimated at p = 5.5 GeV/c
400  Double_t yTrErrMS = 0.00055; // estimated at p = 5.5 GeV/c
401  Double_t xTrErr[nMaxh], yTrErr[nMaxh];
402  for (Int_t i = 0; i < nMaxh; i++) {
403  xTrErr[i] = TMath::Sqrt(xTrErrDet*xTrErrDet+xTrErrMS*xTrErrMS);
404  yTrErr[i] = TMath::Sqrt(yTrErrDet*yTrErrDet+yTrErrMS*yTrErrMS);
405  }
406  Int_t nTr, nCells;
407 
408  // print info on multiple cells
409  //cell = GetCellByGID(3091);
410  //cell->PrintCell("MC");
411  if (kFALSE) {
412  Bool_t recTrack;
413  Int_t nCells, nCells1, nDiffTracks = 0;
414  Int_t nTrackID[10000], TrackID[10000];
415  for (Int_t i = 0; i < 10000; i++) {
416  TrackID[i] = -1;
417  nTrackID[i] = 0;
418  }
419  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
420  layer = GetLayer(iL);
421  nCells = layer->GetNcells();
422 
423  nDiffTracks = 0;
424  for (Int_t i = 0; i < 10000; i++) {
425  TrackID[i] = -1;
426  nTrackID[i] = 0;
427  }
428  nCells1 = 0;
429 
430  for (Int_t iC = 0; iC < nCells; iC++) {
431  cell = layer->GetCell(iC);
432  if (cell->GetTrackGID(0) == cell->GetTrackGID(1)) {
433  nCells1++;
434  recTrack = kTRUE;
435  for (Int_t idt = 0; idt < nDiffTracks; idt++) {
436  if (TrackID[idt] == cell->GetTrackGID(0)) {
437  nTrackID[idt]++;
438  recTrack = kFALSE;
439  break;
440  }
441  }
442  if (recTrack) {
443  TrackID[nDiffTracks] = cell->GetTrackGID(0);
444  nTrackID[nDiffTracks]++;
445  nDiffTracks++;
446  }
447  }
448  }
449 
450  if (nDiffTracks < nCells1) {
451  printf("AnalyzeCells: L %d cells %d , %d diff tracks %d \n",iL,nCells,nCells1,nDiffTracks);
452  }
453  for (Int_t idt = 0; idt < nDiffTracks; idt++) {
454  if (nTrackID[idt] > 1) {
455  for (Int_t iC = 0; iC < nCells; iC++) {
456  cell = layer->GetCell(iC);
457  if (cell->GetTrackGID(0) == TrackID[idt] &&
458  cell->GetTrackGID(1) == TrackID[idt]) {
459  printf("Cell: \n"); cell->PrintCell("MC");
460  }
461  }
462  }
463  }
464 
465  } // end loop layers
466  //return;
467  } // end print info on multiple cells
468 
469  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
470 
471  layer = GetLayer(iL);
472 
473  nCells = layer->GetNcells();
474 
475  for (Int_t iC1 = 0; iC1 < nCells; iC1++) {
476 
477  cell1 = layer->GetCell(iC1);
478  if (cell1->GetStatus() < 1 || cell1->IsMerged()) continue;
479 
480  zc1h1 = cell1->GetHitp1()[2];
481  zc1h2 = cell1->GetHitp2()[2];
482  trkid1h1 = cell1->GetTrackGID(0);
483  trkid1h2 = cell1->GetTrackGID(1);
484 
485  multCell1 = kFALSE;
486  for (Int_t iC2 = (iC1+1); iC2 < nCells; iC2++) {
487 
488  cell2 = layer->GetCell(iC2);
489  if (cell2->GetStatus() < 1 || cell2->IsMerged()) continue;
490 
491  if (cell2->GetLength() != cell1->GetLength()) continue;
492 
493  zc2h1 = cell2->GetHitp1()[2];
494  zc2h2 = cell2->GetHitp2()[2];
495  trkid2h1 = cell2->GetTrackGID(0);
496  trkid2h2 = cell2->GetTrackGID(1);
497 
498  if ((TMath::Abs(zc1h1-zc2h1) > 0.5*fZGap) ||
499  (TMath::Abs(zc1h2-zc2h2) > 0.5*fZGap)) {
500 
501  xc1h1 = cell1->GetHit1()[0];
502  xc1h2 = cell1->GetHit2()[0];
503  xc2h1 = cell2->GetHit1()[0];
504  xc2h2 = cell2->GetHit2()[0];
505  yc1h1 = cell1->GetHit1()[1];
506  yc1h2 = cell1->GetHit2()[1];
507  yc2h1 = cell2->GetHit1()[1];
508  yc2h2 = cell2->GetHit2()[1];
509  /*
510  // angle between the two cells
511  cellAngDif = GetCellAngleDif(cell1,cell2);
512  // intercept point at the mid distance between layers
513  zint = 0.5*(fPlanesZ[cell1->GetLayers()[0]]+fPlanesZ[cell1->GetLayers()[1]]);
514  cellIntDifX = GetCellInterceptX(cell1,zint)-GetCellInterceptX(cell2,zint);
515  cellIntDifY = GetCellInterceptY(cell1,zint)-GetCellInterceptY(cell2,zint);
516  */
517  cellDifX1 = xc1h1-xc2h1;
518  cellDifX2 = xc1h2-xc2h2;
519  cellDifY1 = yc1h1-yc2h1;
520  cellDifY2 = yc1h2-yc2h2;
521 
522  if (trkid1h1 == trkid2h1 && trkid1h2 == trkid2h2) {
523  cellAngDif = GetCellAngleDif(cell1,cell2);
524  hAngDifDup->Fill(cellAngDif);
525  hIntDifXYDup->Fill(cellDifX1,cellDifY1);
526  hIntDifXYDup->Fill(cellDifX2,cellDifY2);
527  //printf("CellIntDif %f %f %f %f \n",cellDifX1,cellDifX2,cellDifY1,cellDifY2);
528  }
529 
530  //if ((cellAngDif < cellAngDifCut) &&
531  // (TMath::Abs(cellIntDifX) < fXCut) &&
532  // (TMath::Abs(cellIntDifY) < fYCut)) {
533  if (TMath::Abs(cellDifX1) < cellIntDifCut &&
534  TMath::Abs(cellDifX2) < cellIntDifCut &&
535  TMath::Abs(cellDifY1) < cellIntDifCut &&
536  TMath::Abs(cellDifY2) < cellIntDifCut) {
537  if (!multCell1) {
538  //printf("Add mult cell1 %d \n",cell1->GetGID());
539  //cell1->PrintCell("MC");
540  multCell1 = kTRUE;
541  multCell->Add(cell1);
542  }
543  //printf("Add mult cell2 %d \n",cell2->GetGID());
544  //cell2->PrintCell("MC");
545  multCell->Add(cell2);
546  }
547  }
548  /*
549  if (trkid1h1 != trkid2h1 && trkid1h2 != trkid2h2) {
550  // angle between the two cells
551  cellAngDif = GetCellAngleDif(cell1,cell2);
552  // intercept point at the mid distance between layers
553  zint = 0.5*(fPlanesZ[cell1->GetLayers()[0]]+fPlanesZ[cell1->GetLayers()[1]]);
554  cellIntDifX = GetCellInterceptX(cell1,zint)-GetCellInterceptX(cell2,zint);
555  cellIntDifY = GetCellInterceptY(cell1,zint)-GetCellInterceptY(cell2,zint);
556  cellIntDif = TMath::Sqrt(cellIntDifX*cellIntDifX+cellIntDifY*cellIntDifY);
557  hAngDifAll->Fill(cellAngDif);
558  hIntDifXYAll->Fill(cellIntDifX,cellIntDifY);
559  }
560  */
561  } // end cell2 loop
562 
563  // merge multi cells
564  if (multCell->GetSize() > 0) {
565 
566  //printf("multCell size %d \n",multCell->GetSize());
567  nTr = 0;
568  zMin = +9999.;
569  zMax = -9999.;
570  //printf("Found %d multi cells.\n",multCell->GetSize());
571  for (Int_t imc = 0; imc < multCell->GetSize(); imc++) {
572  cellM = (AliMFTCACell*)multCell->At(imc);
573  xTr[nTr] = cellM->GetHit1()[0];
574  yTr[nTr] = cellM->GetHit1()[1];
575  zTr[nTr] = cellM->GetHit1()[2];
576  zMin = TMath::Min(zMin,zTr[nTr]);
577  zMax = TMath::Max(zMax,zTr[nTr]);
578  nTr++;
579  xTr[nTr] = cellM->GetHit2()[0];
580  yTr[nTr] = cellM->GetHit2()[1];
581  zTr[nTr] = cellM->GetHit2()[2];
582  zMin = TMath::Min(zMin,zTr[nTr]);
583  zMax = TMath::Max(zMax,zTr[nTr]);
584  nTr++;
585  cellM->SetStatus(0);
586  //printf("Rm cell %d \n",cellM->GetGID());
587  //cellM->PrintCell("MC");
588  }
589  /*
590  printf("Number of points: %d \n",nTr);
591  for (Int_t iTr = 0; iTr < nTr; iTr++) {
592  printf("%d %f %f %f \n",iTr,xTr[iTr],yTr[iTr],zTr[iTr]);
593  }
594  */
595  if (LinFit(nTr,zTr,xTr,xTrErr,ax,axe,bx,bxe) &&
596  LinFit(nTr,zTr,yTr,yTrErr,ay,aye,by,bye)) {
597  hit1[0] = ax*zMin+bx;
598  hit1[1] = ay*zMin+by;
599  hit1[2] = zMin;
600  hit2[0] = ax*zMax+bx;
601  hit2[1] = ay*zMax+by;
602  hit2[2] = zMax;
603  // add a new cell
604  cell = layer->AddCell();
605  cell->SetHits(hit1,hit2,fPlanesZ[cell1->GetLayers()[0]],fPlanesZ[cell1->GetLayers()[1]]);
606  cell->SetLayers(cell1->GetLayers()[0],cell1->GetLayers()[1]);
607  cell->SetStatus(1);
608  cell->SetGID(fCellGID++,trkid1h1,trkid1h2);
609  cell->SetIsMerged();
610  //printf("Add merged cell %d \n",cell->GetGID());
611  //printf("H1: %f %f %f \n",hit1[0],hit1[1],hit1[2]);
612  //printf("H2: %f %f %f \n",hit2[0],hit2[1],hit2[2]);
613  /*
614  for (Int_t imc = 0; imc < multCell->GetSize(); imc++) {
615  cellM = (AliMFTCACell*)multCell->At(imc);
616  cellAngDif = GetCellAngleDif(cell,cellM);
617  hAngDifDup->Fill(cellAngDif);
618  }
619  */
620  } else {
621  printf("No line fit possible!\n");
622  }
623  multCell->Clear();
624  } // end merge multi cells
625 
626  } // end cell1 loop
627 
628  } // end layer loop
629 
630  delete multCell;
631 
632 }
633 
634 //___________________________________________________________________________
636 
637  // create long cells (with one missing layer in the middle) starting from the
638  // short cells which did not find neighbours during the first RunForward()
639 
640  Bool_t prn = kFALSE;
641 
642  AliMFTCALayer *layerL;
643  AliMFTCACell *cellL;
644  Double_t h1[3], h2[3];
645  Int_t iL2, nH2, trackID1, trackID2;
646 
647  Int_t nGapCells = 0;
648 
649  for (Int_t iL1 = 0; iL1 < (fNlayers-1); iL1++) {
650 
651  layerL = GetLayer(iL1);
652 
653  for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) { // Loop over cell in layer iL
654 
655  // short cell
656  cellL = layerL->GetCell(iCL);
657 
658  // attach at right a long cell
659  if ((iL1 < (fNlayers-3)) && (cellL->GetNNbR() == 0)) {
660 
661  h1[0] = cellL->GetHit2()[0];
662  h1[1] = cellL->GetHit2()[1];
663  h1[2] = cellL->GetHit2()[2];
664 
665  iL2 = iL1 + 3;
666  nH2 = GetLayer(iL2)->GetNhits();
667 
668  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
669 
670  h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
671  h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
672  h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
673 
674  TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
675  if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
676 
677  if (!RuleSelect2LayersGap(iL1+1,iL1+2,h1,h2)) continue;
678 
679  if (!RuleSelectCell(h1,h2,iL1+1)) continue;
680 
681  nGapCells++;
682 
683  AliMFTCACell *cell = GetLayer(iL1+1)->AddCell();
684  cell->SetHits(h1,h2,fPlanesZ[iL1+1],fPlanesZ[iL1+3]);
685  cell->SetLayers(iL1+1,iL1+3);
686  cell->SetStatus(1);
687  trackID1 = cellL->GetTrackGID(1);
688  trackID2 = GetLayer(iL1+3)->GetHit(iH2)->GetTrackGID();
689  cell->SetGID(fCellGID++,trackID1,trackID2);
690 
691  if (prn) printf("Create gap (L) cell %d in layer %d. \n",cell->GetGID(),iL1+1);
692 
693  } // end loop hits
694 
695  } // end attach at right a long cell
696 
697  // attach at left a long cell; 1 <> 2
698  if ((iL1 > 1) && (cellL->GetNNbL() == 0)) {
699 
700  h1[0] = cellL->GetHit1()[0];
701  h1[1] = cellL->GetHit1()[1];
702  h1[2] = cellL->GetHit1()[2];
703 
704  iL2 = iL1 - 2;
705  nH2 = GetLayer(iL2)->GetNhits();
706 
707  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
708 
709  h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
710  h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
711  h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
712 
713  TVector3 vec(h1[0]-h2[0],h1[1]-h2[1],h1[2]-h2[2]);
714  if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
715 
716  if (!RuleSelect2LayersGap(iL1-2,iL1-1,h2,h1)) continue;
717 
718  if (!RuleSelectCell(h2,h1,iL1-2)) continue;
719 
720  nGapCells++;
721 
722  AliMFTCACell *cell = GetLayer(iL1-2)->AddCell();
723  cell->SetHits(h2,h1,fPlanesZ[iL1-2],fPlanesZ[iL1]);
724  cell->SetLayers(iL1-2,iL1);
725  cell->SetStatus(1);
726  trackID2 = cellL->GetTrackGID(0);
727  trackID1 = GetLayer(iL1-2)->GetHit(iH2)->GetTrackGID();
728  cell->SetGID(fCellGID++,trackID1,trackID2);
729 
730  if (prn) printf("Create gap (R) cell %d in layer %d. \n",cell->GetGID(),iL1-2);
731 
732  } // end loop hits
733 
734  } // end attach at left a long cell
735 
736  } // end loop cells
737 
738  } // end loop layers
739 
740  printf("Found %d gap cells.\n",nGapCells);
741 
742 }
743 
744 //___________________________________________________________________________
746 
747  // create cells from hits selected in roads
748 
749  Bool_t prn = kFALSE;
750 
751  AliMFTCACell *cell;
752  AliMFTCAHit *hit1, *hit2;
753  Int_t iL1, iL2, nH1, nH2;
754  Int_t iL1min, iL1max;
755  Int_t iL2min, iL2max;
756  Int_t trackID1, trackID2, detElemID1, detElemID2;
757  Bool_t noCell;
758  Double_t h1[3], h2[3], h[3], hx, hy, dR;
759  Int_t mftClsId1, mftClsId2;
760  Int_t nCombi = 0;
761 
762  iL1min = road->GetLayer1();
763  iL1max = road->GetLayer2()-1;
764  //printf("R%d iL1 %d %d \n",ir,iL1min,iL1max);
765 
766  for (iL1 = iL1min; iL1 <= iL1max; iL1++) {
767  //printf("iL1 %d \n",iL1);
768  iL2min = iL1+1;
769  nH1 = road->GetNhitsInLayer(iL1);
770  for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
771  //printf("iH1 %d \n",iH1);
772  hit1 = road->GetHitInLayer(iL1,iH1);
773  iL2max = TMath::Min((iL1+(5-hit1->IsFace())),fNlayers-1);
774  h1[0] = hit1->GetPos()[0];
775  h1[1] = hit1->GetPos()[1];
776  h1[2] = hit1->GetPos()[2];
777  mftClsId1 = hit1->GetMFTClsId();
778  noCell = kTRUE;
779  iL2 = iL2min;
780  while (noCell && (iL2 <= iL2max)) {
781  //printf("iL2 %d \n",iL2);
782  nH2 = road->GetNhitsInLayer(iL2);
783  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
784  hit2 = road->GetHitInLayer(iL2,iH2);
785  h2[0] = hit2->GetPos()[0];
786  h2[1] = hit2->GetPos()[1];
787  h2[2] = hit2->GetPos()[2];
788  mftClsId2 = hit2->GetMFTClsId();
789  nCombi++;
790  if (RuleSelectCell(h1,h2,iL1)) {
791  noCell = kFALSE;
792  cell = GetLayer(iL1)->AddCell();
793  cell->SetHits(h1,h2,fPlanesZ[iL1],fPlanesZ[iL2]);
794  cell->SetMFTClsId(mftClsId1,mftClsId2);
795  cell->SetLayers(iL1,iL2);
796  cell->SetStatus(1);
797  trackID1 = hit1->GetTrackGID();
798  trackID2 = hit2->GetTrackGID();
799  detElemID1 = hit1->GetDetElemID();
800  detElemID2 = hit2->GetDetElemID();
801  cell->SetGID(fCellGID++,trackID1,trackID2);
802  cell->SetDetElemID(detElemID1,detElemID2);
803  road->AddCell(cell);
804  //printf("Cell %5d: L %d-%d H %03d-%03d MC %04d-%04d \n",cell->GetGID(),iL1,iL2,iH1,iH2,trackID1,trackID2);
805  } // end create cell
806  } // end loop iH2
807  iL2++;
808  } // end loop iL2
809  } // end loop iH1
810  } // end loop iL1
811  /*
812  for (Int_t iL = 0; iL < fNlayers; iL++) {
813  printf("%5d %2d %5d \n",ir,iL,road->GetNhitsInLayer(iL));
814  }
815  */
816 
817  if (kFALSE) {
818  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
819  for (Int_t iC = 0; iC < road->GetNcellsInLayer(iL); iC++) {
820  cell = road->GetCellInLayer(iL,iC);
821  printf("L%d,%d-%d CellGID %d MC %d %d \n",iL,cell->GetLayers()[0],cell->GetLayers()[1],cell->GetGID(),cell->GetTrackGID(0),cell->GetTrackGID(1));
822  }
823  }
824  }
825 
826  if (kFALSE) {
827  printf("From %d combinations: \n",nCombi);
828  Long_t nTotCell =0;
829  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
830  printf("Layer %d nr of cells: %d \n",iL,GetLayer(iL)->GetNcells());
831  nTotCell += GetLayer(iL)->GetNcells();
832  }
833  printf("Tot cells: %ld \n",nTotCell);
834 
835  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
836  Int_t nc = road->GetNcellsInLayer(iL);
837  printf("L%1d C%2d \n",iL,nc);
838  }
839 
840  }
841 
842 }
843 
844 //___________________________________________________________________________
845 void AliMFTTrackFinder::CreateCells(Bool_t cv /*Calculate Vertex*/) {
846 
847  Bool_t prn = kFALSE;
848 
849  AliMFTCACell *cell;
850  AliMFTCAHit *hit1, *hit2;
851  Int_t iL1, iL2, nH1, nH2;
852  Int_t iL1min, iL1max;
853  Int_t iL2min, iL2max;
854  Int_t trackID1, trackID2, detElemID1, detElemID2;
855  Bool_t noCell;
856  Double_t h1[3], h2[3], h[3], hx, hy, dR;
857  Int_t nCombi = 0;
858 
859  iL1min = 0;
860  iL1max = fNlayers-2;
861 
862  for (iL1 = iL1min; iL1 <= iL1max; iL1++) {
863  //printf("iL1 %d \n",iL1);
864  iL2min = iL1+1;
865  nH1 = GetLayer(iL1)->GetNhits();
866  for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
867  //printf("iH1 %d \n",iH1);
868  hit1 = GetLayer(iL1)->GetHit(iH1);
869  iL2max = TMath::Min((iL1+(5-hit1->IsFace())),fNlayers-1);
870  h1[0] = hit1->GetPos()[0];
871  h1[1] = hit1->GetPos()[1];
872  h1[2] = hit1->GetPos()[2];
873  noCell = kTRUE;
874  iL2 = iL2min;
875  while (noCell && (iL2 <= iL2max)) {
876  //printf("iL2 %d \n",iL2);
877  nH2 = GetLayer(iL2)->GetNhits();
878  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
879  hit2 = GetLayer(iL2)->GetHit(iH2);
880  h2[0] = hit2->GetPos()[0];
881  h2[1] = hit2->GetPos()[1];
882  h2[2] = hit2->GetPos()[2];
883  nCombi++;
884  if (RuleSelectCell(h1,h2,iL1)) {
885  noCell = kFALSE;
886  cell = GetLayer(iL1)->AddCell();
887  cell->SetHits(h1,h2,fPlanesZ[iL1],fPlanesZ[iL2]);
888  cell->SetLayers(iL1,iL2);
889  cell->SetStatus(1);
890  trackID1 = hit1->GetTrackGID();
891  trackID2 = hit2->GetTrackGID();
892  detElemID1 = hit1->GetDetElemID();
893  detElemID2 = hit2->GetDetElemID();
894  cell->SetGID(fCellGID++,trackID1,trackID2);
895  cell->SetDetElemID(detElemID1,detElemID2);
896  //printf("Cell %5d: L %d-%d H %03d-%03d MC %04d-%04d \n",cell->GetGID(),iL1,iL2,iH1,iH2,trackID1,trackID2);
897  } // end create cell
898  } // end loop iH2
899  iL2++;
900  } // end loop iL2
901  } // end loop iH1
902  } // end loop iL1
903  /*
904  for (Int_t iL = 0; iL < fNlayers; iL++) {
905  printf("%5d %2d %5d \n",ir,iL,road->GetNhitsInLayer(iL));
906  }
907  */
908 
909  if (kFALSE) {
910  printf("From %d combinations: \n",nCombi);
911  Long_t nTotCell =0;
912  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
913  printf("Layer %d nr of cells: %d \n",iL,GetLayer(iL)->GetNcells());
914  nTotCell += GetLayer(iL)->GetNcells();
915  }
916  printf("Tot cells: %ld \n",nTotCell);
917 
918  }
919 
920 }
921 
922 //___________________________________________________________________________
923 void AliMFTTrackFinder::CreateCellsOld(Bool_t cv /*Calculate Vertex*/) {
924 
925  // create only short cells (between two subsequent layers)
926 
927  // print info on multiple cells
928  if (kFALSE) {
929  Bool_t recTrack;
930  Int_t nHits, nDiffTracks = 0;
931  Int_t nTrackID[10000], TrackID[10000];
932  for (Int_t i = 0; i < 10000; i++) {
933  TrackID[i] = -1;
934  nTrackID[i] = 0;
935  }
936  for (Int_t iL = 0; iL < fNlayers; iL++) {
937  nHits = GetLayer(iL)->GetNhits();
938  nDiffTracks = 0;
939  for (Int_t i = 0; i < 10000; i++) {
940  TrackID[i] = -1;
941  nTrackID[i] = 0;
942  }
943  for (Int_t iH = 0; iH < nHits; iH++) {
944  recTrack = kTRUE;
945  for (Int_t idt = 0; idt < nDiffTracks; idt++) {
946  if (TrackID[idt] == GetLayer(iL)->GetHit(iH)->GetTrackGID()) {
947  nTrackID[idt]++;
948  recTrack = kFALSE;
949  break;
950  }
951  }
952  if (recTrack) {
953  TrackID[nDiffTracks] = GetLayer(iL)->GetHit(iH)->GetTrackGID();
954  nTrackID[nDiffTracks]++;
955  nDiffTracks++;
956  }
957  }
958  printf("CreateCells: L %d hits %d diff tracks %d \n",iL,nHits,nDiffTracks);
959  }
960  }
961 
962  Bool_t prn = kFALSE;
963 
964  // loop over the layers
965  Int_t nH1, nH2, trackID1, trackID2, detElemID1, detElemID2;
966  Double_t h1[3], h2[3];
967  Double_t nCombi = 0.;
968 
969  Int_t iL1, iL2;
970 
971  iL1 = 0;
972  while (iL1 < (fNlayers-1)) {
973 
974  iL2 = iL1 + 1;
975 
976  nH1 = GetLayer(iL1)->GetNhits();
977  if (prn) printf("---> 1st Layer L %d with %d hits.\n",iL1,nH1);
978 
979  nH2 = GetLayer(iL2)->GetNhits();
980  if (prn) printf("---> 2nd Layer R %d with %d hits.\n",iL2,nH2);
981 
982  for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
983 
984  h1[0] = GetLayer(iL1)->GetHit(iH1)->GetPos()[0];
985  h1[1] = GetLayer(iL1)->GetHit(iH1)->GetPos()[1];
986  h1[2] = GetLayer(iL1)->GetHit(iH1)->GetPos()[2];
987 
988  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
989 
990  h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
991  h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
992  h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
993 
994  TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
995  if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
996 
997  nCombi += 1.;
998 
999  if (!cv && !RuleSelectCell(h1,h2,iL1)) continue;
1000 
1001  AliMFTCACell *cell = GetLayer(iL1)->AddCell();
1002  cell->SetHits(h1,h2,fPlanesZ[iL1],fPlanesZ[iL2]);
1003  cell->SetLayers(iL1,iL2);
1004  cell->SetStatus(1);
1005  trackID1 = GetLayer(iL1)->GetHit(iH1)->GetTrackGID();
1006  trackID2 = GetLayer(iL2)->GetHit(iH2)->GetTrackGID();
1007  detElemID1 = GetLayer(iL1)->GetHit(iH1)->GetDetElemID();
1008  detElemID2 = GetLayer(iL2)->GetHit(iH2)->GetDetElemID();
1009  cell->SetGID(fCellGID++,trackID1,trackID2);
1010  cell->SetDetElemID(detElemID1,detElemID2);
1011  //cell->PrintCell("FULL");
1012  //printf("Cell nr: %d \n",GetLayer(iL1)->GetNcells());
1013  //cell = GetLayer(iL1)->GetCell(GetLayer(iL1)->GetNcells()-1);
1014  //cell->PrintCell("FULL");
1015 
1016  } // end loop hit in layer 2
1017 
1018  } // end loop hit in layer 1
1019 
1020  if (prn) printf("Create cell %d in layer %d-%d. \n",GetLayer(iL1)->GetNcells(),iL1,GetLayer(iL1)->GetID());
1021 
1022  if (cv) break;
1023 
1024  iL1++;
1025 
1026  } // end loop layer 1
1027 
1028  if (cv) CalculateVertex();
1029 
1030  if (kTRUE || prn) {
1031  printf("From %.0f combinations: \n",nCombi);
1032  Long_t nTotCell =0;
1033  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
1034  printf("Layer %d nr of cells: %d \n",iL,GetLayer(iL)->GetNcells());
1035  nTotCell += GetLayer(iL)->GetNcells();
1036  }
1037  printf("Tot cells: %ld \n",nTotCell);
1038  }
1039 
1040 }
1041 
1042 //___________________________________________________________________________
1044 
1045  hVertZ->Reset();
1046  hVertZa->Reset();
1047 
1048  AliMFTCALayer *layer = GetLayer(0);
1049  AliMFTCACell *cell;
1050  const Double_t *h1, *h2;
1051  Double_t a, b, c, x0, y0, z0;
1052  Double_t n1, n2, n3, n4;
1053  Double_t zmin;
1054 
1055  for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
1056 
1057  cell = layer->GetCell(iC);
1058 
1059  h1 = cell->GetHit1();
1060  h2 = cell->GetHit2();
1061 
1062  x0 = h1[0];
1063  y0 = h1[1];
1064  z0 = h1[2];
1065  a = (h2[0]-h1[0])/(h2[2]-h1[2]);
1066  b = (h2[1]-h1[1])/(h2[2]-h1[2]);
1067  c = 1.;
1068 
1069  n1 = (c*x0-a*z0)/c;
1070  n2 = a/c;
1071  n3 = (c*y0-b*z0)/c;
1072  n4 = b/c;
1073 
1074  zmin = -(n1*n2+n3*n4)/(n2*n2+n4*n4);
1075 
1076  hVertZ->Fill(zmin);
1077  hVertZa->Fill(zmin);
1078 
1079  }
1080 
1081  Float_t zvert = 0., sum = 0., maxsum = 0.;
1082  Int_t bin1, bin2, binW = 2;
1083  Int_t binMin = 1;
1084  Int_t binMax = hVertZ->GetNbinsX();
1085 
1086  hVertZ->Fit("pol1","QW0");
1087  TF1 *f = hVertZ->GetFunction("pol1");
1088 
1089  for (Int_t i = binMin; i <= binMax; i++) {
1090  hVertZ->SetBinContent(i,TMath::Max(0.,hVertZ->GetBinContent(i)-f->Eval(hVertZ->GetBinCenter(i))));
1091  }
1092 
1093  for (Int_t i = binMin; i <= binMax; i++) {
1094  bin1 = TMath::Max(binMin,i-binW);
1095  bin2 = TMath::Min(binMax,i+binW);
1096  sum = hVertZ->Integral(bin1,bin2);
1097  if (sum > maxsum) {
1098  maxsum = sum;
1099  zvert = hVertZ->GetBinCenter(i);
1100  }
1101  //printf("%d %d %d %f %f %f \n",bin1,i,bin2,sum,maxsum,zvert);
1102  }
1103  fZVertCalc = zvert;
1104  printf("Fit vertex z = %f cm \n",fZVertCalc);
1105 
1106  // range = zvert +/- 3 cm
1107  fZVertRange[0] = fZVertCalc - 3.;
1108  fZVertRange[1] = fZVertCalc + 3.;
1109 
1110 }
1111 
1112 //___________________________________________________________________________
1113 void AliMFTTrackFinder::RunForwardR(AliMFTCARoad *road, Int_t& trackGID) {
1114 
1115  Bool_t prn = kFALSE;
1116 
1117  if (prn) AliInfo("Run forward (roads) ==================================== \n");
1118 
1119  AliMFTCALayer *layerL;
1120  AliMFTCALayer *layerR;
1121  AliMFTCACell *cellL;
1122  AliMFTCACell *cellR;
1123  Bool_t stch;
1124  Int_t iL, iR, iter;
1125  Double_t nCombiTot = 0;
1126  Double_t nCombiMatch = 0;
1127  Int_t cellLayers[2];
1128 
1129  fMaxCellStatus = 0;
1130 
1131  iter = 0;
1132 
1133  stch = kTRUE;
1134  while (stch) {
1135 
1136  stch = kFALSE;
1137  iter++;
1138 
1139  for (iL = 0; iL < (fNlayers-2); iL++) { // loop over layers
1140 
1141  for (Int_t iCL = 0; iCL < road->GetNcellsInLayer(iL); iCL++) {
1142 
1143  cellL = road->GetCellInLayer(iL,iCL);
1144  if (cellL->GetStatus() == 0) continue;
1145 
1146  for(Int_t i = 0; i < 2; i++) cellLayers[i] = cellL->GetLayers()[i];
1147 
1148  iR = iL+(cellLayers[1]-cellLayers[0]);
1149  if (iR >= (fNlayers-1))
1150  continue;
1151 
1152  for (Int_t iCR = 0; iCR < road->GetNcellsInLayer(iR); iCR++) {
1153 
1154  cellR = road->GetCellInLayer(iR,iCR);
1155  if (cellR->GetStatus() == 0) continue;
1156 
1157  nCombiTot += 1;
1158 
1159  if ((cellL->GetStatus() == cellR->GetStatus()) &&
1160  RuleSelect(cellL,cellR)) {
1161 
1162  if (prn){
1163  AliInfo(Form("Matching cells: L(%d) cellGID(%d) %d R(%d) cellGID(%d) %d \n",iL,iCL,cellL->GetGID(),iR,iCR,cellR->GetGID()));
1164  }
1165 
1166  nCombiMatch += 1;
1167 
1168  if (iter == 1) {
1169  cellL->AddRightNeighbour(cellR->GetGID());
1170  cellR->AddLeftNeighbour(cellL->GetGID());
1171  }
1172 
1173  cellR->IncrStatus();
1174 
1175  stch = kTRUE;
1176 
1177  } // END : matching cells
1178 
1179  } // END : loop over cells in layer iR
1180 
1181  } // END : loop over cells in layer iL
1182 
1183  } // END : loop over layer iL
1184 
1186 
1187  if (kFALSE || prn) {
1188  AliInfo(Form("Iteration: %5d ----------------- \n",iter));
1189  for (iL = 0; iL < (fNlayers-1); iL++) {
1190  for (Int_t iCL = 0; iCL < road->GetNcellsInLayer(iL); iCL++) {
1191  cellL = road->GetCellInLayer(iL,iCL);
1192  if (cellL->HasNbL() || cellL->HasNbR()) {
1193  printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1194  }
1195  }
1196  }
1197  }
1198 
1199  } // end status change
1200 
1201  if (kFALSE || prn) {
1202  printf("End iteration: ----------------- \n");
1203  for (iL = 0; iL < (fNlayers-1); iL++) {
1204  for (Int_t iCL = 0; iCL < road->GetNcellsInLayer(iL); iCL++) {
1205  cellL = road->GetCellInLayer(iL,iCL);
1206  if (cellL->HasNbL() || cellL->HasNbR()) {
1207  printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1208  }
1209  }
1210  }
1211  }
1212 
1213  if (kFALSE || prn) {
1214 
1215  printf("RunForward after %d iterations, nr of Total combinations %.0f , nr of matched combinations %.0f\n",iter,nCombiTot, nCombiMatch);
1216 
1217  printf("After RunForward max cell status = %d \n",fMaxCellStatus);
1218 
1219  }
1220 
1221  RunBackwardR(road,trackGID);
1222 
1223 }
1224 
1225 //___________________________________________________________________________
1227 
1228  Bool_t prn = kFALSE;
1229 
1230  if (prn) printf("Run forward ==================================== \n");
1231 
1232  AliMFTCALayer *layerL;
1233  AliMFTCALayer *layerR;
1234  AliMFTCACell *cellL;
1235  AliMFTCACell *cellR;
1236  Bool_t stch = kTRUE;
1237  Int_t iL, iR, iter = 0;
1238  Double_t nCombiTot = 0;
1239  Double_t nCombiMatch = 0;
1240 
1241  Double_t nCombiIter, nCombiIterMatch;
1242 
1243  while (stch) {
1244 
1245  nCombiIter = 0;
1246  nCombiIterMatch = 0;
1247 
1248  stch = kFALSE;
1249  iter++;
1250 
1251  for (iL = 0; iL < (fNlayers-2); iL++) { // loop over layers
1252 
1253  layerL = GetLayer(iL);
1254 
1255  if (prn) printf("L %d cells %d \n",iL,layerL->GetNcells());
1256 
1257  for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) { // Loop over cell in layer iL
1258 
1259  cellL = layerL->GetCell(iCL);
1260  if (cellL->GetStatus() == 0) continue;
1261 
1262  Int_t cellLayers[2];
1263  for(Int_t i = 0; i < 2; i++) cellLayers[i] = cellL->GetLayers()[i];
1264 
1265  iR = iL+(cellLayers[1] - cellLayers[0]);
1266  if (iR < (fNlayers-1))
1267  layerR = GetLayer(iR);
1268  else
1269  continue;
1270 
1271  // vertex selection here ?
1272  //if (!RuleSelectCell(cellL)) continue;
1273 
1274  for (Int_t iCR = 0; iCR < layerR->GetNcells(); iCR++) { // Loop over cell in layer iL
1275 
1276  cellR = layerR->GetCell(iCR);
1277  if (cellR->GetStatus() == 0) continue;
1278 
1279  // vertex selection here ?
1280  //if (!RuleSelectCell(cellR)) continue;
1281 
1282  nCombiTot += 1;
1283  nCombiIter += 1;
1284 
1285  if ((cellL->GetStatus() == cellR->GetStatus()) && RuleSelect(cellL,cellR)) { // If Cells are matching in angles and have equal status values
1286 
1287  if (prn){
1288  printf("Matching cells: L cellGID %d R cellGID %d \n",cellL->GetGID(),cellR->GetGID());
1289  printf("Layer L %d cell %d - Layer R %d cell %d \n",iL,iCL,iR,iCR);
1290  }
1291  nCombiMatch += 1;
1292  nCombiIterMatch += 1;
1293 
1294  if (iter == 1) {
1295  cellL->AddRightNeighbour(cellR->GetGID());
1296  cellR->AddLeftNeighbour(cellL->GetGID());
1297  }
1298 
1299  cellR->IncrStatus();
1300 
1301  stch = kTRUE;
1302 
1303  } // END : matching cells
1304 
1305  } // END : loop over cell in layer iL-1
1306 
1307  } // END : loop over cell layer iL
1308  } // END : loop over layers
1309 
1310  if (prn) {
1311  printf("Iteration: %5d nr of combinations %.0f match %.0f \n",iter,nCombiIter,nCombiIterMatch);
1312  }
1313 
1314  UpdateCellStatus();
1315 
1316  if (prn) {
1317  printf("Iteration: %5d ----------------- \n",iter);
1318  for (iL = 0; iL < (fNlayers-1); iL++) {
1319  layerL = GetLayer(iL);
1320  for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) {
1321  cellL = layerL->GetCell(iCL);
1322  if (cellL->HasNbL() || cellL->HasNbR()) {
1323  printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1324  }
1325  }
1326  }
1327  }
1328 
1329  } // end status change
1330 
1331  if (prn) {
1332  printf("End iteration: ----------------- \n");
1333  for (iL = 0; iL < (fNlayers-1); iL++) {
1334  layerL = GetLayer(iL);
1335  for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) {
1336  cellL = layerL->GetCell(iCL);
1337  if (cellL->HasNbL() || cellL->HasNbR()) {
1338  printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1339  }
1340  }
1341  }
1342  }
1343 
1344  if (kTRUE || prn) {
1345 
1346  printf("RunForward after %d iterations, nr of Total combinations %.0f , nr of matched combinations %.0f\n",iter,nCombiTot, nCombiMatch);
1347 
1348  printf("After RunForward max cell status = %d \n",fMaxCellStatus);
1349 
1350  }
1351 
1352 }
1353 
1354 //___________________________________________________________________________
1355 void AliMFTTrackFinder::RunBackwardR(AliMFTCARoad *road, Int_t& trackGID) {
1356 
1357  Bool_t prn = kFALSE;
1358 
1359  if (prn) printf("Run backward road %d ==================================== \n",road->GetID());
1360 
1361  if (fMaxCellStatus == 1) return; // only isolated cells
1362 
1363  Double_t chisqNbLprev, cellAngDif, cellAngDifPrev;
1364  Int_t addCellIdToTrack, iNbLchiSqMin, iCellAngDifMin, nHitSta;
1365 
1366  AliMFTCALayer *layerR;
1367  AliMFTCACell *cellR;
1368  AliMFTCACell *cellL;
1369  AliMFTCACell *cell;
1371 
1372  Bool_t addCellToTrack, hitSta[5];
1373 
1374  Int_t minStartLayer = 6;
1375  Int_t maxStartLayer = 8;
1376 
1377  for (Int_t startLayer = maxStartLayer; startLayer >= minStartLayer; startLayer--) {
1378 
1379  if (prn) printf("Start layer %d \n",startLayer);
1380 
1381  for (Int_t iCR = 0; iCR < road->GetNcellsInLayer(startLayer); iCR++) {
1382 
1383  cellR = road->GetCellInLayer(startLayer,iCR);
1384 
1385  if (cellR->GetStatus() == 0) continue;
1386  if (cellR->IsUsed()) continue;
1387  if (cellR->GetStatus() < (fMinTrackLength-1)) continue;
1388 
1389  if (prn) printf("Create new track %d \n",trackGID);
1390 
1391  track = AddTrack(trackGID++);
1392  track->SetStartLayer(startLayer);
1393  track->AddCell(cellR);
1394  cellR->SetUsed(kTRUE);
1395 
1396  // add cells to new track
1397  addCellToTrack = kTRUE;
1398  while (addCellToTrack) {
1399 
1400  cellR = road->GetCellByGID(track->GetLastCellGID()); // !!!
1401 
1402  addCellToTrack = kFALSE;
1403 
1404  // find the left neighbor giving the smalles chisquare
1405  iNbLchiSqMin = 0;
1406  chisqNbLprev = 0.;
1407 
1408  // find the left neighbor giving the smallest deviation
1409  cellAngDifPrev = -1.;
1410  iCellAngDifMin = 0;
1411 
1412  // do a first loop to check all possible associations
1413  for (Int_t iNNbL = 0; iNNbL < cellR->GetNNbL(); iNNbL++) {
1414 
1415  cellL = road->GetCellByGID(cellR->GetNbLgid(iNNbL));
1416 
1417  if (kFALSE && prn) {
1418  printf("To track %d attach cell GID {L}: %d - TrackID of this cell : %d - %d, Status %d \n",track->GetGID(),cellL->GetGID(),cellL->GetTrackGID(0),cellL->GetTrackGID(1),cellL->GetStatus());
1419  }
1420 
1421  if (cellL->GetStatus() == 0) continue;
1422  if (cellL->IsUsed()) continue;
1423  if (cellL->GetStatus() != (cellR->GetStatus()-1)) continue;
1424 
1425  // ... smallest deviation
1426  cellAngDif = GetCellAngleDif(cellL,cellR);
1427  if (cellAngDifPrev < 0.) {
1428  cellAngDifPrev = cellAngDif;
1429  } else {
1430  if (cellAngDif < cellAngDifPrev) {
1431  iCellAngDifMin = iNNbL;
1432  }
1433  }
1434 
1435  // ... smallest chisquare
1436  if (track->AddCellToChiSq(cellL) < chisqNbLprev) {
1437  iNbLchiSqMin = iNNbL;
1438  }
1439 
1440  chisqNbLprev = track->AddCellToChiSq(cellL);
1441 
1442  } // END : left neighbour loop
1443 
1444  //if (cellR->GetNNbL() > 1) {
1445  // printf("%d %d \n",iNbLchiSqMin,iCellAngDifMin);
1446  //}
1447 
1448  // use the angular deviation instead of the chisquare
1449  //iNbLchiSqMin = iCellAngDifMin;
1450 
1451  // do a second loop and take the good association of cells
1452  if (cellR->GetNNbL() > 0) {
1453  cellL = road->GetCellByGID(cellR->GetNbLgid(iNbLchiSqMin));
1454  addCellToTrack = kTRUE;
1455  addCellIdToTrack = cellR->GetNbLgid(iNbLchiSqMin);
1456  cellL = road->GetCellByGID(addCellIdToTrack);
1457  track->AddCell(cellL);
1458  cellL->SetUsed(kTRUE);
1459  if (prn) {
1460  printf("To track %d attach cell GID {L}: %d - TrackID of this cell : %d - %d, Status %d \n",track->GetGID(),cellL->GetGID(),cellL->GetTrackGID(0),cellL->GetTrackGID(1),cellL->GetStatus());
1461  }
1462  }
1463 
1464  } // END : addCellToTrack
1465 
1466  // check again track length
1467  for (Int_t j = 0; j < 5; j++) hitSta[j] = kFALSE;
1468  for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1469  cell = GetCellByGID(track->GetCellGID(iCell));
1470  hitSta[cell->GetLayers()[0]/2] = kTRUE;
1471  hitSta[cell->GetLayers()[1]/2] = kTRUE;
1472  }
1473  nHitSta = 0;
1474  for (Int_t i = 0; i < 5; i++) {
1475  if (hitSta[i]) nHitSta++;
1476  }
1477  if (nHitSta < fMinTrackLength) {
1478  for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1479  cell = track->GetCell(iCell);
1480  cell->SetUsed(kFALSE);
1481  }
1482  RemoveLastTrack();
1483  trackGID--;
1484  } else {
1485  road->SetGood();
1486  }
1487 
1488  } // END : startLayer cells loop
1489 
1490  } // END : startLayer loop
1491 
1492 }
1493 
1494 //___________________________________________________________________________
1496 
1497  Bool_t prn = kFALSE;
1498 
1499  if (prn) printf("Run backward ==================================== \n");
1500 
1501  if (fMaxCellStatus == 1) return; // only isolated cells
1502 
1503  Double_t chisqNbLprev, cellAngDif, cellAngDifPrev;
1504  Int_t addCellIdToTrack, iNbLchiSqMin, iCellAngDifMin, trackGID = 0;
1505 
1506  AliMFTCALayer *layerR;
1507  AliMFTCACell *cellR;
1508  AliMFTCACell *cellL;
1509  AliMFTCACell *cell;
1511 
1512  Bool_t addCellToTrack;
1513 
1514  Int_t minStartLayer = fMinTrackLength-2;
1515  Int_t maxStartLayer = (fNlayers-1)-1;
1516 
1517  for (Int_t startLayer = maxStartLayer; startLayer >= minStartLayer; startLayer--) {
1518 
1519  if (prn) printf("Start layer %d \n",startLayer);
1520 
1521  layerR = GetLayer(startLayer);
1522 
1523  for (Int_t iCR = 0; iCR < layerR->GetNcells(); iCR++) {
1524  cellR = layerR->GetCell(iCR);
1525  if (cellR->GetStatus() == 0) continue;
1526  if (cellR->IsUsed()) continue;
1527  if (cellR->GetStatus() >= (fMinTrackLength-1)) {
1528  if (prn) printf("Create new track %d \n",trackGID);
1529  track = AddTrack(trackGID++);
1530  track->SetStartLayer(startLayer);
1531  //track->AddCellGID(cellR->GetGID());
1532  track->AddCell(cellR);
1533  cellR->SetUsed(kTRUE);
1534  if (prn) {
1535  printf("To track %d attach cell GID: %d - TrackID of this cell : %d - %d\n",track->GetGID(),cellR->GetGID(),cellR->GetTrackGID(0),cellR->GetTrackGID(1));
1536  }
1537  // add cells to new track
1538  addCellToTrack = kTRUE;
1539  while (addCellToTrack) {
1540  cellR = GetCellByGID(track->GetLastCellGID());
1541  addCellToTrack = kFALSE;
1542  iNbLchiSqMin = 0;
1543  chisqNbLprev = 0.;
1544  cellAngDifPrev = -1.;
1545  iCellAngDifMin = 0;
1546  for (Int_t iNNbL = 0; iNNbL < cellR->GetNNbL(); iNNbL++) {
1547  cellL = GetCellByGID(cellR->GetNbLgid(iNNbL));
1548  if (cellL->GetStatus() == 0) continue;
1549  if (cellL->IsUsed()) continue;
1550  if (cellL->GetStatus() == (cellR->GetStatus()-1)) {
1551  cellAngDif = GetCellAngleDif(cellL,cellR);
1552  if (cellAngDifPrev < 0.) {
1553  cellAngDifPrev = cellAngDif;
1554  } else {
1555  if (cellAngDif < cellAngDifPrev) {
1556  iCellAngDifMin = iNNbL;
1557  }
1558  }
1559  if (track->AddCellToChiSq(cellL) < chisqNbLprev) {
1560  iNbLchiSqMin = iNNbL;
1561  }
1562  }
1563  chisqNbLprev = track->AddCellToChiSq(cellL);
1564  } // END : left neighbour loop
1565  //if (cellR->GetNNbL() > 1) {
1566  // printf("%d %d \n",iNbLchiSqMin,iCellAngDifMin);
1567  //}
1568  // use the angular deviation instead of the chisquare
1569  //iNbLchiSqMin = iCellAngDifMin;
1570  if (cellR->GetNNbL() > 0) {
1571  cellL = GetCellByGID(cellR->GetNbLgid(iNbLchiSqMin));
1572  addCellToTrack = kTRUE;
1573  addCellIdToTrack = cellR->GetNbLgid(iNbLchiSqMin);
1574  //track->AddCellGID(addCellIdToTrack);
1575  cellL = GetCellByGID(addCellIdToTrack);
1576  track->AddCell(cellL);
1577  cellL->SetUsed(kTRUE);
1578  if (prn) {
1579  printf("To track %d attach cell GID: %d - TrackID of this cell : %d - %d\n",track->GetGID(),cellL->GetGID(),cellL->GetTrackGID(0),cellL->GetTrackGID(1));
1580  }
1581  }
1582  } // END : addCellToTrack
1583 
1584  // check again track length
1585  if (track->GetNcells() < (fMinTrackLength-1)) {
1586  for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1587  cell = track->GetCell(iCell);
1588  cell->SetUsed(kFALSE);
1589  }
1590  RemoveLastTrack();
1591  trackGID--;
1592  }
1593 
1594  } // END : create new track
1595 
1596  } // END : startLayer cells loop
1597 
1598  } // END : startLayer loop
1599 
1600 }
1601 
1602 //___________________________________________________________________________
1604  AliCodeTimerAuto("",0);
1605 
1606  Bool_t prn = kFALSE;
1607 
1608  AliInfo(Form("Filtering %d tracks",GetNtracks()));
1609 
1610  Int_t nTrkC = 0, nTrkG = 0, nTrkF = 0, nTrkN = 0;
1611 
1613  AliMFTCACell *cell, *celln;
1614  Int_t ndof, nptr, cellGID, cellGIDn, cellGIDprev = -1, nTrkSplitEnd = 0;
1615  const Int_t nMaxh = 100;
1616  Double_t xTr[nMaxh], yTr[nMaxh], zTr[nMaxh];
1617  Double_t a, ae, b, be, x0, xS, y0, yS, zmin, chisqx, chisqy;
1618  Double_t trkPhi, trkThe;
1619  Bool_t splitTrack, recTrack, cleanTrack, goodTrack, fakeTrack, noiseTrack;
1620  const Int_t nTrackMax = 10000;
1621  Int_t nrHitTrackID[nTrackMax], idHitTrackID[nTrackMax], nDiffTracks;
1622  Int_t idMaxHitsTrack, nMaxHits, nNoisyPix;
1623  for (Int_t i = 0; i < nTrackMax; i++) {
1624  nrHitTrackID[i] = 0;
1625  idHitTrackID[i] = -2;
1626  }
1627  Int_t nrAliMFTCATrackID[nTrackMax], idAliMFTCATrackID[nTrackMax], nDiffAliMFTCATracks = 0;
1628  for (Int_t i = 0; i < nTrackMax; i++) {
1629  nrAliMFTCATrackID[i] = 0;
1630  idAliMFTCATrackID[i] = -2;
1631  }
1632  Double_t xTrErrDet = 0.0028/TMath::Sqrt(12.);
1633  Double_t yTrErrDet = 0.0028/TMath::Sqrt(12.);
1634  Double_t xTrErrMS = 0.00055; // estimated at p = 5.5 GeV/c
1635  Double_t yTrErrMS = 0.00055; // estimated at p = 5.5 GeV/c
1636  Double_t xTrErr[nMaxh], yTrErr[nMaxh];
1637  for (Int_t i = 0; i < nMaxh; i++) {
1638  xTrErr[i] = TMath::Sqrt(xTrErrDet*xTrErrDet+xTrErrMS*xTrErrMS);
1639  yTrErr[i] = TMath::Sqrt(yTrErrDet*yTrErrDet+yTrErrMS*yTrErrMS);
1640  }
1641  fErrX = xTrErr[0];
1642  fErrY = yTrErr[0];
1643 
1644  Int_t nTotalHits = 0, nCleanTotalHits = 0;
1645 
1646  for (Int_t iTrack = 0; iTrack < GetNtracks(); iTrack++) {
1647  track = GetTrack(iTrack);
1648  nDiffTracks = 0;
1649  nptr = 0;
1650  for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1651  cellGID = track->GetCellGID(iCell);
1652  cell = GetCellByGID(cellGID);
1653  //printf("Track %3d Cell %5d \n",iTrack,cellGID);
1654  // tracks split from the first (highest status) cell
1655  if (iCell == 0) {
1656  if (cellGIDprev >= 0) {
1657  if (cellGID == cellGIDprev) {
1658  // this is a split track
1659  nTrkSplitEnd++;
1660  splitTrack = kTRUE;
1661  } else {
1662  splitTrack = kFALSE;
1663  }
1664  }
1665  cellGIDprev = cellGID;
1666  }
1667  if (prn) {
1668  printf("Cell %4d from track %d ",cellGID,track->GetGID());
1669  printf("(%4d %4d %d) \n",cell->GetTrackGID(0),cell->GetTrackGID(1),cell->GetGID());
1670  printf("2: %f %f %f \n",cell->GetHit2()[0],cell->GetHit2()[1],cell->GetHit2()[2]);
1671  printf("1: %f %f %f \n",cell->GetHit1()[0],cell->GetHit1()[1],cell->GetHit1()[2]);
1672  }
1673  // extract hit x,y,z
1674  if (nptr == 0) {
1675  xTr[nptr] = cell->GetHit2()[0];
1676  yTr[nptr] = cell->GetHit2()[1];
1677  zTr[nptr] = cell->GetHit2()[2];
1678  nptr++;
1679  xTr[nptr] = cell->GetHit1()[0];
1680  yTr[nptr] = cell->GetHit1()[1];
1681  zTr[nptr] = cell->GetHit1()[2];
1682  nptr++;
1683  } else {
1684  xTr[nptr] = cell->GetHit1()[0];
1685  yTr[nptr] = cell->GetHit1()[1];
1686  zTr[nptr] = cell->GetHit1()[2];
1687  nptr++;
1688  }
1689  // count the generated tracks which contribute to this reconstructed track
1690  for (Int_t ihc = 0; ihc < 2; ihc++) {
1691  recTrack = kTRUE;
1692  for (Int_t idt = 0; idt < nDiffTracks; idt++) {
1693  if (idHitTrackID[idt] == cell->GetTrackGID(ihc)) {
1694  nrHitTrackID[idt]++;
1695  recTrack = kFALSE;
1696  break;
1697  }
1698  }
1699  if (recTrack) {
1700  idHitTrackID[nDiffTracks] = cell->GetTrackGID(ihc);
1701  nrHitTrackID[nDiffTracks] = 1;
1702  nDiffTracks++;
1703  }
1704  } // cell hits
1705  /*
1706  // debug
1707  if (kFALSE || pTot[cell->GetTrackGID(0)] > 4.0) {
1708  RuleSelectCell(cell);
1709  if (iCell < (track->GetNcells()-1)) {
1710  cellGIDn = track->GetCellGID(iCell+1);
1711  celln = GetCellByGID(cellGIDn);
1712  SetDebug(1);
1713  RuleSelect(celln,cell);
1714  SetDebug(0);
1715  }
1716  }
1717  //
1718  */
1719  } // end cell loop
1720  // assert quality of the track
1721  cleanTrack = goodTrack = fakeTrack = noiseTrack = kFALSE;
1722  if (nDiffTracks == 1 && idHitTrackID[0] >= 0) {
1723  cleanTrack = kTRUE;
1724  idMaxHitsTrack = idHitTrackID[0];
1725  } else {
1726  nNoisyPix = 0;
1727  nMaxHits = 0;
1728  for (Int_t idt = 0; idt < nDiffTracks; idt++) {
1729  if (idHitTrackID[idt] == -1) {
1730  nNoisyPix = nrHitTrackID[idt];
1731  } else if (nMaxHits < nrHitTrackID[idt]) {
1732  nMaxHits = nrHitTrackID[idt];
1733  idMaxHitsTrack = idHitTrackID[idt];
1734  }
1735  }
1736  if (GetNDet() == 5) {
1737  // allow one fake hit
1738  if (nMaxHits >= (2*track->GetNcells())-1) {
1739  goodTrack = kTRUE;
1740  } else {
1741  if (nNoisyPix > 0) noiseTrack = kTRUE;
1742  else fakeTrack = kTRUE;
1743  }
1744  }
1745  if (GetNDet() == 5*2) {
1746  // allow two fake hits
1747  if (nMaxHits >= (2*track->GetNcells())-2) {
1748  goodTrack = kTRUE;
1749  } else {
1750  if (nNoisyPix > 0) noiseTrack = kTRUE;
1751  else fakeTrack = kTRUE;
1752  }
1753  }
1754  if (GetNDet() == 6) {
1755  // allow one fake hit
1756  if (nMaxHits >= (2*track->GetNcells())-1) {
1757  goodTrack = kTRUE;
1758  } else {
1759  if (nNoisyPix > 0) noiseTrack = kTRUE;
1760  else fakeTrack = kTRUE;
1761  }
1762  }
1763  } // end assert track quality
1764  nTotalHits += nptr;
1765  if (cleanTrack) {
1766  // count the duplicated clean tracks
1767  recTrack = kTRUE;
1768  for (Int_t idt = 0; idt < nDiffAliMFTCATracks; idt++) {
1769  if (idAliMFTCATrackID[idt] == idMaxHitsTrack) {
1770  nrAliMFTCATrackID[idt]++;
1771  recTrack = kFALSE;
1772  break;
1773  }
1774  }
1775  if (recTrack) {
1776  idAliMFTCATrackID[nDiffAliMFTCATracks] = idMaxHitsTrack;
1777  nrAliMFTCATrackID[nDiffAliMFTCATracks] = 1;
1778  nDiffAliMFTCATracks++;
1779  }
1780  track->SetMCflag(1);
1781  track->SetMCindex(idMaxHitsTrack);
1782  nTrkC++;
1783  nCleanTotalHits += nptr;
1784  }
1785  if (goodTrack) {
1786  track->SetMCflag(2);
1787  nTrkG++;
1788  }
1789  if (fakeTrack) {
1790  track->SetMCflag(3);
1791  nTrkF++;
1792  }
1793  if (noiseTrack) {
1794  track->SetMCflag(4);
1795  nTrkN++;
1796  }
1797  // linear regression
1798  //printf("Fit line with %d points.\n",nptr);
1799  if (LinFit(nptr,zTr,xTr,xTrErr,a,ae,b,be)) {
1800  x0 = b; xS = a;
1801  if (LinFit(nptr,zTr,yTr,yTrErr,a,ae,b,be)) {
1802  y0 = b; yS = a;
1803  chisqx = 0.;
1804  chisqy = 0.;
1805  for (Int_t iptr = 0; iptr < nptr; iptr++) {
1806  //printf("%d %f %f %f \n",iptr,xTr[iptr],yTr[iptr],zTr[iptr]);
1807  chisqx += (xTr[iptr]-(xS*zTr[iptr]+x0))*(xTr[iptr]-(xS*zTr[iptr]+x0))/(xTrErr[iptr]*xTrErr[iptr]);
1808  chisqy += (yTr[iptr]-(yS*zTr[iptr]+y0))*(yTr[iptr]-(yS*zTr[iptr]+y0))/(yTrErr[iptr]*yTrErr[iptr]);
1809  }
1810  // track phi and theta
1811  trkPhi = trkThe = 0.;
1812  if (TMath::Abs(xS) > 0.) {
1813  trkPhi = TMath::ATan(yS/xS);
1814  // put the correct signs
1815  if (xS < 0. && yS > 0.) {
1816  trkPhi = TMath::Pi()+trkPhi;
1817  }
1818  if (xS < 0. && yS < 0.) {
1819  trkPhi = TMath::Pi()+trkPhi;
1820  }
1821  if (xS > 0. && yS < 0.) {
1822  trkPhi = TMath::TwoPi()+trkPhi;
1823  }
1824  if (TMath::Abs(TMath::Sin(trkPhi)) > 0.) {
1825  trkThe = TMath::ATan(yS/TMath::Sin(trkPhi));
1826  // put the correct signs
1827  trkThe = TMath::Abs(trkThe); // is always smaller than 90deg
1828  trkPhi *= TMath::RadToDeg();
1829  trkThe *= TMath::RadToDeg();
1830  }
1831  }
1832  // calculate DCA with the beam axis
1833  zmin = -(x0*xS+y0*yS)/(xS*xS+yS*yS);
1834  track->SetTheta(trkThe);
1835  track->SetPhi(trkPhi);
1836  if (fCalcVertex) {
1837  track->SetVertX(x0+xS*fZVertCalc);
1838  track->SetVertY(y0+yS*fZVertCalc);
1839  track->SetVertZ(fZVertCalc);
1840  } else {
1841  track->SetVertX(x0+xS*fZVertDet);
1842  track->SetVertY(y0+yS*fZVertDet);
1843  track->SetVertZ(fZVertDet);
1844  }
1845  ndof = nptr-2;
1846  track->SetChiSqX(chisqx/(Double_t)ndof);
1847  track->SetChiSqY(chisqy/(Double_t)ndof);
1848  } // yz fit
1849  } // xz fit
1850  //printf("End fit.\n");
1851  } // end track loop
1852 
1853  fNDifTracks = nDiffAliMFTCATracks;
1854 
1855  AliInfo(Form("Track found -> C: %5d G: %5d F: %5d N: %5d Dif: %5d TotalHits: %d Clean: %d",nTrkC,nTrkG,nTrkF,nTrkN,nDiffAliMFTCATracks,nTotalHits,nCleanTotalHits));
1856 
1857 }
1858 
1859 //___________________________________________________________________________
1860 void AliMFTTrackFinder::DrawTracks(Double_t *pTot, Double_t *Theta) {
1861 
1862  Bool_t prn = kFALSE;
1863 
1865  AliMFTCACell *cell;
1866  Int_t cellGID, trackGID, nTrackS = 0;
1867  Bool_t single[10000], hitFromNoisyPix, hitFromDiffTrack;
1868  Int_t nGoodCell, firstHitTrackID, idTrack[10000], nGoodTracks = 0;
1869  for (Int_t i = 0; i < 10000; i++) idTrack[i] = -1;
1870 
1871  printf("Draw %d tracks. \n",GetNtracks());
1872  /*
1873  Double_t x[nDet], y[nDet], z[nDet], a, ae, b, be;
1874  Double_t errx[nDet], erry[nDet];
1875  for (Int_t i = 0; i < nDet; i++) {
1876  errx[i] = det[i]->GetSigmaX();
1877  erry[i] = det[i]->GetSigmaY();
1878  }
1879  */
1880  //Double_t chisq;
1881  for (Int_t iT = 0; iT < GetNtracks(); iT++) {
1882  // cell content in a track
1883  single[iT] = kTRUE;
1884  track = GetTrack(iT);
1885  cellGID = track->GetCellGID(0);
1886  cell = GetCellByGID(cellGID);
1887  trackGID = cell->GetTrackGID(0);
1888  single[iT] &= (cell->GetTrackGID(1) == trackGID);
1889  nGoodCell = 0;
1890  for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1891  cellGID = track->GetCellGID(iC);
1892  cell = GetCellByGID(cellGID);
1893  single[iT] &= (cell->GetTrackGID(0) == trackGID);
1894  single[iT] &= (cell->GetTrackGID(1) == trackGID);
1895  if (single[iT]) nGoodCell++;
1896  //x[iC] = cell->GetHit1()[0];
1897  //y[iC] = cell->GetHit1()[1];
1898  //z[iC] = cell->GetHit1()[2];
1899  }
1900  //____________________________________________________________________
1901  /*
1902  chisq = 0.;
1903  if (LinFit(track->GetNcells(),z,x,errx,a,ae,b,be)) {
1904  for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1905  chisq += (x[iC]-(a*z[iC]+b))*(x[iC]-(a*z[iC]+b));
1906  }
1907  #ifdef HOUGH
1908  hRTheXZ->Fill(90.-TMath::ATan(a)*TMath::RadToDeg(),b*TMath::Sin(TMath::ATan(a)));
1909  #endif
1910  }
1911  if (LinFit(track->GetNcells(),z,y,erry,a,ae,b,be)) {
1912  for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1913  chisq += (y[iC]-(a*z[iC]+b))*(y[iC]-(a*z[iC]+b));
1914  }
1915  #ifdef HOUGH
1916  hRTheYZ->Fill(90.-TMath::ATan(a)*TMath::RadToDeg(),b*TMath::Sin(TMath::ATan(a)));
1917  #endif
1918  }
1919  if (prn) printf("ChiSq = %e \n",chisq);
1920  */
1921  //____________________________________________________________________
1922  if (fDebug > 0) {
1923  hNGoodCell->Fill(nGoodCell);
1924  }
1925  if (single[iT]) {
1926  nTrackS++;
1927  } else {
1928  if (prn || fDebug > 0) {
1929  for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1930  cellGID = track->GetCellGID(iC);
1931  cell = GetCellByGID(cellGID);
1932  if (prn) printf("Track %4d Cell %2d TrackGID %5d %5d \n",iT,iC,cell->GetTrackGID(0),cell->GetTrackGID(1));
1933  }
1934  }
1935  }
1936  if (prn) printf("Track %d with %d cells: \n",iT,track->GetNcells());
1937  hitFromNoisyPix = kFALSE;
1938  hitFromDiffTrack = kFALSE;
1939  firstHitTrackID = -1;
1940  for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1941  cellGID = track->GetCellGID(iC);
1942  cell = GetCellByGID(cellGID);
1943  if (iC == 0) firstHitTrackID = cell->GetTrackGID(0);
1944  if (prn) printf("\t%5d from track: ",cellGID);
1945  if (prn) printf("%5d %5d \n",cell->GetTrackGID(0),cell->GetTrackGID(1));
1946  if (cell->GetTrackGID(0) == -1 || cell->GetTrackGID(1) == -1)
1947  hitFromNoisyPix = kTRUE;
1948  else if (cell->GetTrackGID(0) != firstHitTrackID ||
1949  cell->GetTrackGID(1) != firstHitTrackID)
1950  hitFromDiffTrack = kTRUE;
1951 
1952  //cell->PrintCell("FULL");
1953  }
1954  if (hitFromDiffTrack) {
1955  hTrackType->Fill(1);
1956  //PrintTrack(iT);
1957  } else if (hitFromNoisyPix) {
1958  hTrackType->Fill(2);
1959  //PrintTrack(iT);
1960  } else {
1961  if (idTrack[firstHitTrackID] >= 0) {
1962  hTrackType->Fill(3);
1963  printf("Double track %6d from %6d and %6d \n",firstHitTrackID,iT,idTrack[firstHitTrackID]);
1964  PrintTrack(iT);
1965  PrintTrack(idTrack[firstHitTrackID]);
1966  } else {
1967  idTrack[firstHitTrackID] = iT;
1968  hTrackType->Fill(0);
1969  if (track->GetNcells() >= (fMinTrackLength-1)) nGoodTracks++;
1970  //PrintTrack(iT);
1971  //printf("%d \n",track->GetNcells());
1972  }
1973  }
1974  } // end loop tracks
1975 
1976  printf("... %d single. \n",nTrackS);
1977  printf("... %d (%d) good \n",(Int_t)hTrackType->GetBinContent(1),nGoodTracks);
1978 
1979 }
1980 
1981 //___________________________________________________________________________
1983 
1984  AliMFTCALayer *layer;
1985  AliMFTCACell *cell;
1986 
1987  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
1988  layer = GetLayer(iL);
1989  for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
1990  cell = layer->GetCell(iC);
1991  if (gid == cell->GetGID()) return cell;
1992  }
1993  }
1994 
1995  return 0;
1996 
1997 }
1998 
1999 //___________________________________________________________________________
2000 Bool_t AliMFTTrackFinder::RuleSelect(AliMFTCACell *cellL, AliMFTCACell *cellR) { // Are cells matching in angle ?
2001  /*
2002  if (0) {
2003  // ideal
2004  if (cellL->GetTrackGID(0) != cellL->GetTrackGID(1)) return kFALSE;
2005  if (cellR->GetTrackGID(0) != cellR->GetTrackGID(1)) return kFALSE;
2006  if (cellL->GetTrackGID(1) != cellR->GetTrackGID(0)) return kFALSE;
2007  }
2008  */
2009  //if (cellL->GetLayers()[1] != cellR->GetLayers()[0]) {
2010  // printf("Neighbor cells do not touch the same common layer!\n");
2011  //}
2012  Int_t layer = cellL->GetLayers()[1];
2013 
2014  TVector3 *segL_ = cellL->GetSeg();
2015  TVector3 *segR_ = cellR->GetSeg();
2016 
2017  TVector3 segL = TVector3(*segL_);
2018  TVector3 segR = TVector3(*segR_);
2019 
2020  const Double_t *hitL[2];
2021  const Double_t *hitR[2];
2022 
2023  hitL[0] = cellL->GetHit1();
2024  hitL[1] = cellL->GetHit2();
2025  hitR[0] = cellR->GetHit1();
2026  hitR[1] = cellR->GetHit2();
2027 
2028  //printf("%f %f %f - %f %f %f \n",hitL[1][0],hitL[1][1],hitL[1][2],hitR[0][0],hitR[0][1],hitR[0][2]);
2029 
2030  Double_t dx, dy, a;
2031  dx = (hitL[1][0]-hitR[0][0])*(hitL[1][0]-hitR[0][0]); // Distance in X direction between the 2 hits of the two different cells in the same layer
2032  dy = (hitL[1][1]-hitR[0][1])*(hitL[1][1]-hitR[0][1]); // Distance in Y direction between the 2 hits of the two different cells in the same layer
2033  dx = (dx > 0.) ? (TMath::Sqrt(dx)) : 0.;
2034  dy = (dy > 0.) ? (TMath::Sqrt(dy)) : 0.;
2035  a = (segL.Angle(segR))*TMath::RadToDeg(); // Angle between the segments of each cell
2036 
2037  if (fDebug > 0) {
2038  hDA[cellL->GetLayers()[1]]->Fill(a);
2039  hDXY[cellL->GetLayers()[1]]->Fill(dx*1.E4,dy*1.E4); // in microns
2040  //hDXY[cellL->GetLayers()[1]]->Fill(dx,dy);
2041  /*
2042  printf("--------------------\n");
2043  segL.Print();
2044  segR.Print();
2045  TVector3 segL1 = TVector3(hitL[1][0]-hitL[0][0],
2046  hitL[1][1]-hitL[0][1],
2047  hitL[1][2]-hitL[0][2]);
2048  TVector3 segR1 = TVector3(hitR[1][0]-hitR[0][0],
2049  hitR[1][1]-hitR[0][1],
2050  hitR[1][2]-hitR[0][2]);
2051  segL1.Print();
2052  segR1.Print();
2053  printf("%f %f %f \n",hitL[0][0],hitL[0][1],hitL[0][2]);
2054  printf("%f %f %f - %f %f %f \n",hitL[1][0],hitL[1][1],hitL[1][2],hitR[0][0],hitR[0][1],hitR[0][2]);
2055  printf("%f %f %f \n",hitR[1][0],hitR[1][1],hitR[1][2]);
2056  */
2057  /*
2058  if (cellL->GetTrackGID(0) == cellL->GetTrackGID(1) &&
2059  cellL->GetTrackGID(1) == cellR->GetTrackGID(0) &&
2060  cellR->GetTrackGID(0) == cellR->GetTrackGID(1)) {
2061  hDA[cellL->GetLayers()[1]]->Fill(a);
2062  hDXY[cellL->GetLayers()[1]]->Fill(dx,dy);
2063  }
2064  */
2065  }
2066  /*
2067  if (a < fACutN[layer]) {
2068  printf("dx, dy, a : %f %f %f \n",dx,dy,a);
2069  }
2070  */
2071  //printf("RS: %f %f %f %f %f %f \n",dx,fXCut,dy,fYCut,a,fACutN[layer]);
2072  if ((dx > fXCut) || (dy > fYCut) || (a > fACutN[layer])) return kFALSE;
2073 
2074  return kTRUE;
2075 
2076 }
2077 
2078 //___________________________________________________________________________
2079 Bool_t AliMFTTrackFinder::RuleSelect2LayersGap(Int_t iL1, Int_t iL2, Double_t *hit1, Double_t *hit2) { // Are Cell formed by hit1 and hit2 compatible with any cell of layers 1 and 2 ?
2080 
2081  Bool_t prn = kFALSE;
2082 
2083  Bool_t findCompatibleCells = kFALSE;
2084 
2085  AliMFTCALayer *layer1 = 0;
2086  if (iL1 >= 0) layer1 = GetLayer(iL1);
2087  else return kFALSE;
2088 
2089  AliMFTCALayer *layer2 = 0;
2090  if (iL2 >= 0) layer2 = GetLayer(iL2);
2091 
2092  for (Int_t iCell = 0; iCell < layer1->GetNcells(); iCell++) {
2093 
2094  AliMFTCACell * cell = layer1->GetCell(iCell);
2095  const Double_t *hit = cell->GetHit1();
2096 
2097  // Distance in X direction between the 2 hits of the two different cells in the same layer
2098  Double_t dx = (hit[0]-hit1[0])*(hit[0]-hit1[0]);
2099 
2100  // Distance in Y direction between the 2 hits of the two different cells in the same layer
2101  Double_t dy = (hit[1]-hit1[1])*(hit[1]-hit1[1]);
2102 
2103  //Double_t radius = (dx+dy > 0.) ? (TMath::Sqrt(dx+dy)) : 0.;
2104 
2105  if ((TMath::Abs(dx) > fXCut) || (TMath::Abs(dy) > fYCut)) continue;
2106 
2107  TVector3 *seg1_ = cell->GetSeg();
2108  TVector3 seg1 = TVector3(*seg1_);
2109 
2110  TVector3 seg2(hit2[0]-hit1[0],hit2[1]-hit1[1],hit2[2]-hit1[2]);
2111  Double_t a = (seg1.Angle(seg2))*TMath::RadToDeg(); // Angle between the segments of each cell
2112 
2113  //printf(Form("Angle = %f\n",a));
2114  if (fDebug > 0) hAngleCells->Fill(a);
2115 
2116  if (prn) {
2117  printf("Compare with cell %d in layer %d \n",iCell,layer1->GetID());
2118  }
2119 
2120  if (a < 0.1) findCompatibleCells = kTRUE;
2121 
2122  }
2123 
2124  if (iL2 < 0) return (!findCompatibleCells);
2125 
2126  for (Int_t iCell = 0; iCell < layer2->GetNcells(); iCell++) {
2127 
2128  AliMFTCACell * cell = layer2->GetCell(iCell);
2129  const Double_t *hit = cell->GetHit2();
2130 
2131  // Distance in X direction between the 2 hits of the two different cells in the same layer
2132  Double_t dx = (hit[0]-hit2[0])*(hit[0]-hit2[0]);
2133 
2134  // Distance in Y direction between the 2 hits of the two different cells in the same layer
2135  Double_t dy = (hit[1]-hit2[1])*(hit[1]-hit2[1]);
2136 
2137  //Double_t radius = (dx+dy > 0.) ? (TMath::Sqrt(dx+dy)) : 0.;
2138 
2139  if ((TMath::Abs(dx) > fXCut) || (TMath::Abs(dy) > fYCut)) continue;
2140 
2141  TVector3 *seg1_ = cell->GetSeg();
2142  TVector3 seg1 = TVector3(*seg1_);
2143 
2144  TVector3 seg2(hit2[0]-hit1[0],hit2[1]-hit1[1],hit2[2]-hit1[2]);
2145  Double_t a = (seg1.Angle(seg2))*TMath::RadToDeg(); // Angle between the segments of each cell
2146 
2147  //printf(Form("Angle = %f\n",a));
2148  if (fDebug > 0) hAngleCells->Fill(a);
2149 
2150  if (prn) {
2151  printf("Compare with cell %d in layer %d \n",iCell,layer2->GetID());
2152  }
2153 
2154  if (a < 0.1) findCompatibleCells = kTRUE;
2155 
2156  }
2157 
2158  return (!findCompatibleCells);
2159 
2160 }
2161 
2162 
2163 
2164 //___________________________________________________________________________
2165 Bool_t AliMFTTrackFinder::RuleSelectCell(AliMFTCACell *cell) { // Look if segment pointing to the vertex
2166 
2167  Int_t layer = cell->GetLayers()[0];
2168 
2169  TVector3 *seg_ = cell->GetSeg();
2170 
2171  TVector3 seg = TVector3(*seg_);
2172 
2173  const Double_t *hit;
2174 
2175  hit = cell->GetHit1();
2176 
2177  Double_t vert[3] = { 0., 0., 0. };
2178  if (fCalcVertex) vert[2] = fZVertCalc;
2179  else vert[2] = fZVertDet;
2180 
2181  Double_t a;
2182  TVector3 segV;
2183 
2184  segV = TVector3(hit[0]-vert[0],hit[1]-vert[1],hit[2]-vert[2]);
2185  a = (seg.Angle(segV))*TMath::RadToDeg();
2186 
2187  hThetaCells->Fill(seg.Theta()*TMath::RadToDeg());
2188  hDAv[layer]->Fill(a);
2189 
2190  if ( a > fACutV[layer]) return kFALSE;
2191 
2192  return kTRUE;
2193 
2194 }
2195 
2196 //___________________________________________________________________________
2197 Bool_t AliMFTTrackFinder::RuleSelectCell(Double_t *h1, Double_t *h2, Int_t iL1, TF1 *f, Bool_t acalc) {
2198 
2199  // Look if segment pointing to the vertex (using directly the hits)
2200 
2201  Double_t a, av[2], acut;
2202  TVector3 segV, segVdet;
2203 
2204  TVector3 seg = TVector3(h2[0]-h1[0], h2[1]-h1[1], h2[2]-h1[2]);
2205 
2206  Double_t vert[3] = { 0., 0., 0. };
2207  if (fCalcVertex) vert[2] = fZVertCalc;
2208  else vert[2] = fZVertDet;
2209 
2210  segVdet = TVector3(h1[0]-vert[0],h1[1]-vert[1],h1[2]-vert[2]);
2211  a = (seg.Angle(segVdet))*TMath::RadToDeg();
2212 
2213  //hDAv[iL1]->Fill(a);
2214 
2215  if (acalc) {
2216  acut = f->Eval(180.-segVdet.Theta()*TMath::RadToDeg());
2217  } else {
2218  acut = fACutV[iL1];
2219  }
2220  //printf("%f %f %d \n",180.-segVdet.Theta()*TMath::RadToDeg(),acut,acalc);
2221  //printf("RSC: %f %f \n",a,acut);
2222 
2223  if ( a > acut) {
2224  AliDebug(2,"Cell NOT Selected");
2225  return kFALSE;
2226  }
2227  AliDebug(2,"Cell Selected");
2228 
2229  return kTRUE;
2230 
2231 }
2232 
2233 //___________________________________________________________________________
2235 
2236  new ((*fTracks)[fNtracks++]) AliMFTCATrack();
2237  AliMFTCATrack *track = (AliMFTCATrack*)fTracks->At(fTracks->GetLast());
2238  track->SetGID(gid);
2239 
2240  return track;
2241 
2242 }
2243 
2244 //___________________________________________________________________________
2246 
2247  // create new track and copy
2248 
2249  new ((*fTracks)[fNtracks++]) AliMFTCATrack();
2250  AliMFTCATrack *track = (AliMFTCATrack*)fTracks->At(fTracks->GetLast());
2251  track->SetGID(gid);
2252 
2253  track->SetStartLayer(trk.GetStartLayer());
2254 
2255  for (Int_t iC = 0; iC < trk.GetNcells(); iC++) {
2256  //track->AddCellGID(trk.GetCellGID(iC));
2257  track->AddCell(trk.GetCell(iC));
2258  }
2259 
2260  return track;
2261 
2262 }
2263 
2264 //___________________________________________________________________________
2266 
2267  AliMFTCARoad *road;
2268  AliMFTCACell *cell;
2269 
2270  for (Int_t ir = 0; ir < fNRoads; ir++) {
2271  road = GetRoad(ir);
2272  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
2273  for (Int_t iC = 0; iC < road->GetNcellsInLayer(iL); iC++) {
2274  cell = road->GetCellInLayer(iL,iC);
2275  cell->UpdateStatus();
2276  fMaxCellStatus = TMath::Max(fMaxCellStatus,cell->GetStatus());
2277  }
2278  }
2279  }
2280 
2281 }
2282 
2283 //___________________________________________________________________________
2285 
2286  AliMFTCALayer *layer;
2287  AliMFTCACell *cell;
2288 
2289  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
2290  layer = GetLayer(iL);
2291  for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
2292  cell = layer->GetCell(iC);
2293  cell->UpdateStatus();
2294  fMaxCellStatus = TMath::Max(fMaxCellStatus,cell->GetStatus());
2295  }
2296  }
2297 
2298 }
2299 
2300 //___________________________________________________________________________
2302 
2303  AliMFTCATrack *track = GetTrack(id);
2304  Int_t cellGID;
2305  AliMFTCACell *cell;
2306  printf("Track:\t%6d \n",id);
2307  for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
2308  cellGID = track->GetCellGID(iC);
2309  cell = GetCellByGID(cellGID);
2310  printf("cell\t%6d gid\t%6d \tfrom track: ",iC,cellGID);
2311  printf("\t%6d\t%6d \tin layers: ",cell->GetTrackGID(0),cell->GetTrackGID(1));
2312  printf("\t%1d\t%1d\n",cell->GetLayers()[0],cell->GetLayers()[1]);
2313  }
2314 
2315 }
2316 
2317 //___________________________________________________________________________
2319 
2320  AliMFTCALayer *layer;
2321  AliMFTCACell *cell;
2322 
2323  for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
2324  layer = GetLayer(iL);
2325  printf("LayerID %d \n",layer->GetID());
2326  for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
2327  cell = layer->GetCell(iC);
2328  cell->PrintCell("FULL");
2329  }
2330  }
2331 
2332 }
2333 
2334 //___________________________________________________________________________
2336 
2337 // TCanvas *ca1 = new TCanvas("CA1","",50,50,800,800);
2338 // TCanvas *ca2 = new TCanvas("CA2","",50,50,800,800);
2339 // TCanvas *ca3 = new TCanvas("CA3","",50,50,800,800);
2340 // ca1->Divide(3,2);
2341 // ca2->Divide(3,2);
2342 // for (Int_t i = 0; i<fNlayers; i++) {
2343 // ca1->cd(i+1);
2344 // hDA[i]->Draw();
2345 // ca2->cd(i+1);
2346 // hDXY[i]->Draw("colz");
2347 // }
2348 // ca3->Divide(1,2);
2349 // ca3->cd(1);
2350 // hNGoodCell->Draw();
2351 // ca3->cd(2);
2352 // hTrackType->Draw();
2353 
2354 }
2355 
2356 //___________________________________________________________________________
2358 
2359  printf("==== AliMFTTrackFinder::PrintParam ====\n");
2360 
2361  printf("Number of layers: %d \n",fNlayers);
2362  printf("Use the TrackFinder: %d \n",fUseTF);
2363  printf("Layer thickness in X0: %5.3f \n",fThick);
2364  printf("Pixel noise: %4.2e \n",fPixelNoise);
2365  printf("Add noise: %d \n",fAddNoise);
2366  printf("fMinTrackLength: %d \n",fMinTrackLength);
2367  printf("fXCut: %6.4f cm\n",fXCut);
2368  printf("fYCut: %6.4f cm\n",fYCut);
2369  printf("fMaxSegAngle: %4.1f deg\n",fMaxSegAngle);
2370  for (Int_t i = 0; i < fNlayers; i++) {
2371  printf("fACutV[%d]: %4.2f \n",i,fACutV[i]);
2372  }
2373  for (Int_t i = 0; i < fNlayers; i++) {
2374  printf("fACutN[%d]: %4.2f \n",i,fACutN[i]);
2375  }
2376  for (Int_t i = 0; i < fNlayers; i++) {
2377  printf("PlaneDetEff[%d]: %4.2f \n",i,fPlaneDetEff[i]);
2378  }
2379  printf("Calculate vertex: %d \n",fCalcVertex);
2380 
2381  printf("CPU time: %f seconds\n",fCPUTime);
2382  printf("Real time: %f seconds\n",fRealTime);
2383 
2384  printf("============================\n");
2385 
2386 }
2387 
2388 //___________________________________________________________________________
2390 
2391  new ((*fRoads)[fNRoads++]) AliMFTCARoad();
2392  AliMFTCARoad *road = (AliMFTCARoad*)fRoads->At(fRoads->GetLast());
2393 
2394  return road;
2395 
2396 }
2397 
2398 //___________________________________________________________________________
2400  AliCodeTimerAuto("",0);
2401 
2402  Bool_t prn = kFALSE;
2403 
2404  // planes: 0, 1, 2, ..., 9 (layer ... becomes plane)
2405  // rules for combining first/last plane in a road:
2406  // 0 with 6, 7, 8, 9
2407  // 1 with 6, 7, 8, 9
2408  // 2 with 8, 9
2409  // 3 with 8, 9
2410 
2411  Int_t iPla1Min = 0, iPla1Max = 3;
2412  Int_t iPla2Min[4] = { 6, 6, 6, 6 };
2413  Int_t iPla2Max[4] = { 9, 9, 7, 7 };
2414 
2415  Int_t nH1, nH2, nH;
2416  Int_t roadLen, nRoads = 0, trackGID = 0;
2417  Double_t h1[3], h2[3], h[3], hx, hy, dR;
2418 
2419  Double_t dRmin = 0.0400;
2420 
2421  AliMFTCAHit *hit1, *hit2, *hit, *htmp;
2422  AliMFTCARoad *road, *road1;
2423  Bool_t hitSta[5] = { kFALSE, kFALSE, kFALSE, kFALSE, kFALSE };
2424  Bool_t isUsed, newRoad;
2425 
2426  for (Int_t iL1 = iPla1Min; iL1 <= iPla1Max; iL1++) {
2427 
2428  for (Int_t iL2 = iPla2Max[iL1]; iL2 >= iPla2Min[iL1]; iL2--) {
2429 
2430  // see AliMFTCARoad::SetLength
2431  roadLen = iL2/2-iL1/2;
2432 
2433  nH1 = GetLayer(iL1)->GetNhits();
2434  if(prn) AliInfo(Form("nH1 = %d ",nH1));
2435  for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
2436  //AliInfo(Form("iH1 = %d ",iH1));
2437 
2438  hit1 = GetLayer(iL1)->GetHit(iH1);
2439 
2440  if (hit1->IsUsed()) continue;
2441  /*
2442  // check if it belongs to a good longer road previously found
2443  isUsed = kFALSE;
2444  for (Int_t i = 0; i < hit1->GetNRoads(); i++) {
2445  road1 = GetRoad(hit1->GetInRoad(i));
2446  if (road1->IsGood()) {
2447  if (road1->GetLength() >= roadLen) {
2448  isUsed = kTRUE;
2449  break;
2450  }
2451  }
2452  }
2453  //if (isUsed) continue;
2454  */
2455  if (prn)
2456  printf("Hit1: %d %d %d %d \n",hit1->GetLayer(),hit1->GetID(),hit1->GetDetElemID(),hit1->GetTrackGID());
2457 
2458  h1[0] = hit1->GetPos()[0];
2459  h1[1] = hit1->GetPos()[1];
2460  h1[2] = hit1->GetPos()[2];
2461 
2462  nH2 = GetLayer(iL2)->GetNhits();
2463  if(prn) AliInfo(Form("nH2 = %d ",nH2));
2464 
2465  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
2466 
2467  hit2 = GetLayer(iL2)->GetHit(iH2);
2468 
2469  if (hit2->IsUsed()) continue;
2470  /*
2471  // check if it belongs to a good longer road previously found
2472  isUsed = kFALSE;
2473  for (Int_t i = 0; i < hit2->GetNRoads(); i++) {
2474  road1 = GetRoad(hit2->GetInRoad(i));
2475  if (road1->IsGood()) {
2476  if (road1->GetLength() >= roadLen) {
2477  isUsed = kTRUE;
2478  break;
2479  }
2480  }
2481  }
2482  //if (isUsed) continue;
2483  */
2484  if (prn)
2485  printf("Hit2: %d %d %d %d \n",hit2->GetLayer(),hit2->GetID(),hit2->GetDetElemID(),hit2->GetTrackGID());
2486 
2487  h2[0] = hit2->GetPos()[0];
2488  h2[1] = hit2->GetPos()[1];
2489  h2[2] = hit2->GetPos()[2];
2490 
2491  TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
2492  if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
2493 
2494  if (!RuleSelectCell(h1,h2,iL1)) continue;
2495 
2496  newRoad = kTRUE;
2497  for (Int_t i = 0; i < 5; i++) hitSta[i] = kFALSE;
2498 
2499  for (Int_t iL = (iL1+1); iL <= (iL2-1); iL++) {
2500 
2501  nH = GetLayer(iL)->GetNhits();
2502 
2503  for (Int_t iH = 0; iH < nH; iH++) {
2504 
2505  hit = GetLayer(iL)->GetHit(iH);
2506 
2507  if (hit->IsUsed()) continue;
2508 
2509  h[0] = hit->GetPos()[0];
2510  h[1] = hit->GetPos()[1];
2511  h[2] = hit->GetPos()[2];
2512 
2513  hx = h1[0] + (h2[0]-h1[0])*(h[2]-h1[2])/(h2[2]-h1[2]);
2514  hy = h1[1] + (h2[1]-h1[1])*(h[2]-h1[2])/(h2[2]-h1[2]);
2515 
2516  dR = TMath::Sqrt((hx-h[0])*(hx-h[0])+(hy-h[1])*(hy-h[1]));
2517  if (dR >= dRmin) continue;
2518  /*
2519  // check if it belongs to a good longer road previously found
2520  isUsed = kFALSE;
2521  for (Int_t i = 0; i < hit->GetNRoads(); i++) {
2522  AliMFTCARoad *road1 = GetRoad(hit->GetInRoad(i));
2523  if (road1->IsGood()) {
2524  if (road1->GetLength() > roadLen) {
2525  isUsed = kTRUE;
2526  break;
2527  }
2528  }
2529  }
2530  //if (isUsed) continue;
2531  */
2532  if (prn)
2533  printf("Hit: %d %d %d %d \n",hit->GetLayer(),hit->GetID(),hit->GetDetElemID(),hit->GetTrackGID());
2534 
2535  if (newRoad) {
2536  //AliInfo("Adding new road");
2537  road = AddRoad();
2538  road->SetID(nRoads++);
2539  road->AddHit(hit1);
2540  road->AddHit(hit2);
2541 
2542  hit1->SetInRoad(road->GetID());
2543  hit2->SetInRoad(road->GetID());
2544 
2545  hitSta[iL1/2] = kTRUE;
2546  hitSta[iL2/2] = kTRUE;
2547 
2548  road->SetLength(iL1,iL2);
2549 
2550  newRoad = kFALSE;
2551 
2552  }
2553 
2554  road->AddHit(hit);
2555  hit->SetInRoad(road->GetID());
2556  hitSta[iL/2] = kTRUE;
2557 
2558  } // end loop hits middle plane
2559 
2560  } // end loop middle plane
2561 
2562  // count the number of hit stations (disks)
2563  if (!newRoad) {
2564  Int_t nHitSta = 0;
2565  for (Int_t i = 0; i < 5; i++)
2566  if (hitSta[i]) nHitSta++;
2567  road->SetNHitSta(nHitSta);
2568  if (nHitSta >= fMinTrackLength) {
2569  CreateCellsR(road);
2570  RunForwardR(road,trackGID);
2571  if (road->IsGood()) {
2572  for (Int_t j = 0; j < fNlayers; j++) {
2573  for (Int_t k = 0; k < road->GetNhitsInLayer(j); k++) {
2574  hit = road->GetHitInLayer(j,k);
2575  //printf("%d %d %d \n",hit->GetLayer(),hit->GetID(),GetLayer(hit->GetLayer())->GetNhits());
2576  htmp = GetLayer(hit->GetLayer())->GetHit(hit->GetID());
2577  htmp->SetUsed();
2578  if (prn)
2579  printf("Hit used: %d %d %d %d \n",hit->GetLayer(),hit->GetID(),hit->GetDetElemID(),hit->GetTrackGID());
2580  }
2581  }
2582  }
2583  }
2584  }
2585 
2586  } // end loop hits last plane
2587 
2588  } // end loop hits first plane
2589 
2590  } // end loop last plane
2591 
2592  } // end loop first plane
2593  /*
2594  Int_t nRoadsGood = 0, nTotalHits = 0;
2595  for (Int_t i = 0; i < nRoads; i++) {
2596  road = GetRoad(i);
2597  if (road->IsGood()) {
2598  //printf("Road %d with %d hits.\n",i,road->GetNhits());
2599  nRoadsGood++;
2600  nTotalHits += road->GetNhits();
2601  }
2602  }
2603  printf("Found %d roads %d good with %d hits. \n",nRoads,nRoadsGood,nTotalHits);
2604  */
2605 }
2606 
2607 //___________________________________________________________________________
2609  AliCodeTimerAuto("",0);
2610 
2611  Bool_t prn = kFALSE;
2612 
2613  Double_t h1[3], h2[3], h[3], hx, hy;
2614  Double_t htr1[3], htr2[3];
2615  const Int_t nMaxh = 100;
2616  Double_t trX[nMaxh], trY[nMaxh], trZ[nMaxh];
2617  Int_t lay[nMaxh], trkid[nMaxh], hit[nMaxh], detid[nMaxh];
2618  Int_t nH1, nH2, nH, iL1, iL2, lay1, lay2, trkid1, trkid2, nptr;
2619  Int_t hit1, hit2, detid1, detid2, nHitSta;
2620  Int_t mftClsId1, mftClsId2;
2621  Int_t trackGID = 0;
2622  AliMFTCACell *cell;
2624 
2625  Double_t dR, dRmin, dRcut = 0.0100;
2626 
2627  TF1 *fACutF = new TF1("fACutF","[0]+x*[1]",0.,1.);
2628  Float_t cut170 = 0.6; // cut [deg] at theta = 170 deg
2629  Float_t cut177 = 0.3; // cut [deg] at theta = 177 deg
2630  Float_t par[2];
2631  par[1] = (cut177-cut170)/(177.-170.);
2632  par[0] = cut170-par[1]*170.;
2633  fACutF->SetParameter(0,par[0]);
2634  fACutF->SetParameter(1,par[1]);
2635 
2636  Bool_t addHit, selCAgapCell;
2637  Bool_t hitSta[5];
2638  Bool_t seed = kTRUE;
2639 
2640  Int_t step = 0;
2641 
2642  iL1 = 0;
2643 
2644  while (seed) {
2645 
2646  if (step == 0) {
2647  iL2 = fNlayers-1;
2648  } else {
2649  iL2--;
2650  }
2651 
2652  step++;
2653 
2654  if (iL2 < iL1 + (fMinTrackLength-1)) {
2655  iL1++;
2656  if (iL1 > fNlayers - (fMinTrackLength-1)) break;
2657  step = 0;
2658  continue;
2659  }
2660 
2661  if (prn) printf("iL1,iL2 %d %d \n",iL1,iL2);
2662  //continue;
2663 
2664  nH1 = GetLayer(iL1)->GetNhits();
2665  nH2 = GetLayer(iL2)->GetNhits();
2666  if (prn) printf("nH1 %d nH2 %d\n",nH1,nH2);
2667  for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
2668 
2669  if (GetLayer(iL1)->GetHit(iH1)->IsUsed()) continue;
2670 
2671 
2672  h1[0] = GetLayer(iL1)->GetHit(iH1)->GetPos()[0];
2673  h1[1] = GetLayer(iL1)->GetHit(iH1)->GetPos()[1];
2674  h1[2] = GetLayer(iL1)->GetHit(iH1)->GetPos()[2];
2675 
2676  if (prn) printf("H1 %d (%.1f,%.1f,%.1f)\n",iH1,h1[0],h1[1],h1[2]);
2677 
2678  for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
2679 
2680  if (GetLayer(iL2)->GetHit(iH2)->IsUsed()) continue;
2681 
2682 
2683  h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
2684  h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
2685  h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
2686 
2687  TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
2688 
2689  if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
2690  //if (prn) printf("Theta =%f (max = %f)\n",vec.Theta(),fMaxSegAngle*TMath::DegToRad());
2691 
2692  //if (!RuleSelectCell(h1,h2,iL1,fACutF,kTRUE)) continue;
2693  if (!RuleSelectCell(h1,h2,iL1)) continue;
2694  if (prn) printf("H2 %d (%.1f,%.1f,%.1f)\n",iH2,h2[0],h2[1],h2[2]);
2695 
2696  // this is a seed connecting the first and the last layers
2697 
2698  for (Int_t i = 0; i < 5; i++) hitSta[i] = kFALSE;
2699  nHitSta = 0;
2700 
2701  hitSta[iL1/2] = kTRUE;
2702  hitSta[iL2/2] = kTRUE;
2703 
2704  nptr = 0;
2705 
2706  trX[nptr] = h1[0];
2707  trY[nptr] = h1[1];
2708  trZ[nptr] = h1[2];
2709  lay[nptr] = iL1;
2710  hit[nptr] = iH1;
2711  trkid[nptr] = GetLayer(iL1)->GetHit(iH1)->GetTrackGID();
2712  detid[nptr] = GetLayer(iL1)->GetHit(iH1)->GetDetElemID();
2713  nptr++;
2714 
2715  for (Int_t iL = (iL1+1); iL <= (iL2-1); iL++) {
2716 
2717  nH = GetLayer(iL)->GetNhits();
2718 
2719  if (prn) printf("L %d nH %d \n",iL,nH);
2720 
2721  dRmin = dRcut;
2722  addHit = kFALSE;
2723 
2724  for (Int_t iH = 0; iH < nH; iH++) {
2725 
2726  //if (prn) printf("H %d \n",iH);
2727 
2728  if (GetLayer(iL)->GetHit(iH)->IsUsed()) continue;
2729 
2730  h[0] = GetLayer(iL)->GetHit(iH)->GetPos()[0];
2731  h[1] = GetLayer(iL)->GetHit(iH)->GetPos()[1];
2732  h[2] = GetLayer(iL)->GetHit(iH)->GetPos()[2];
2733 
2734  hx = h1[0] + (h2[0]-h1[0])*(h[2]-h1[2])/(h2[2]-h1[2]);
2735  hy = h1[1] + (h2[1]-h1[1])*(h[2]-h1[2])/(h2[2]-h1[2]);
2736 
2737  dR = TMath::Sqrt((hx-h[0])*(hx-h[0])+(hy-h[1])*(hy-h[1]));
2738  if (dR >= dRmin) continue;
2739  AliDebug(1,Form("Hit %d added",iH));
2740  hitSta[iL/2] = kTRUE;
2741 
2742  dRmin = dR;
2743 
2744  trX[nptr] = h[0];
2745  trY[nptr] = h[1];
2746  trZ[nptr] = h[2];
2747  lay[nptr] = iL;
2748  hit[nptr] = iH;
2749  trkid[nptr] = GetLayer(iL)->GetHit(iH)->GetTrackGID();
2750  detid[nptr] = GetLayer(iL)->GetHit(iH)->GetDetElemID();
2751 
2752  addHit = kTRUE;
2753 
2754  } // loop hits intermediate layer
2755 
2756  if (addHit) nptr++;
2757 
2758  } // loop intermediate layers
2759 
2760  trX[nptr] = h2[0];
2761  trY[nptr] = h2[1];
2762  trZ[nptr] = h2[2];
2763  lay[nptr] = iL2;
2764  hit[nptr] = iH2;
2765  trkid[nptr] = GetLayer(iL2)->GetHit(iH2)->GetTrackGID();
2766  detid[nptr] = GetLayer(iL2)->GetHit(iH2)->GetDetElemID();
2767  nptr++;
2768 
2769  // create cells and tracks
2770  //
2771  // calculate how many stations have hit(s)
2772  AliDebug(2,Form("nptr %d fMinTrackLength %d",nptr,fMinTrackLength));
2773  if (nptr < fMinTrackLength) continue;
2774  nHitSta = 0;
2775  for (Int_t i = 0; i < 5; i++) {
2776  if (hitSta[i]) nHitSta++;
2777  }
2778  AliDebug(2,Form("nHitSta %d fMinTrackLength %d",nHitSta,fMinTrackLength));
2779  if (nHitSta < fMinTrackLength) continue;
2780  /*
2781  // as in CA impose gap cells with no more than one layer missing
2782  selCAgapCell = kTRUE;
2783  for (Int_t i = 0; i < (nptr-1); i++) {
2784  if (lay[i+1]-lay[i] > 2) {
2785  selCAgapCell = kFALSE;
2786  break;
2787  }
2788  }
2789  if (!selCAgapCell) continue;
2790  */
2791  // create track
2792  AliDebug(1,Form("Adding Track %d ",trackGID));
2793  track = AddTrack(trackGID++);
2794  track->SetStartLayer(iL1);
2795 
2796  // loop over hits in reverse order, like in RunBackward
2797  for (Int_t iptr = (nptr-1); iptr >= 1; iptr--) {
2798 
2799  trkid1 = trkid[iptr-1];
2800  lay1 = lay[iptr-1];
2801  hit1 = hit[iptr-1];
2802  detid1 = detid[iptr-1];
2803 
2804  htr1[0] = GetLayer(lay1)->GetHit(hit1)->GetPos()[0];
2805  htr1[1] = GetLayer(lay1)->GetHit(hit1)->GetPos()[1];
2806  htr1[2] = GetLayer(lay1)->GetHit(hit1)->GetPos()[2];
2807  mftClsId1 = GetLayer(lay1)->GetHit(hit1)->GetMFTClsId();
2808 
2809  GetLayer(lay1)->GetHit(hit1)->SetUsed();
2810 
2811  trkid2 = trkid[iptr];
2812  lay2 = lay[iptr];
2813  hit2 = hit[iptr];
2814  detid2 = detid[iptr];
2815 
2816  htr2[0] = GetLayer(lay2)->GetHit(hit2)->GetPos()[0];
2817  htr2[1] = GetLayer(lay2)->GetHit(hit2)->GetPos()[1];
2818  htr2[2] = GetLayer(lay2)->GetHit(hit2)->GetPos()[2];
2819  mftClsId2 = GetLayer(lay2)->GetHit(hit2)->GetMFTClsId();
2820 
2821  GetLayer(lay2)->GetHit(hit2)->SetUsed();
2822 
2823  // create cells
2824  cell = GetLayer(lay[iptr-1])->AddCell();
2825  cell->SetHits(htr1,htr2,fPlanesZ[lay1],fPlanesZ[lay2]);
2826  cell->SetMFTClsId(mftClsId1,mftClsId2);
2827  cell->SetLayers(lay1,lay2);
2828  cell->SetStatus(1);
2829  cell->SetGID(fCellGID++,trkid1,trkid2);
2830  cell->SetDetElemID(detid1,detid2);
2831 
2832  // add cell to track
2833  track->AddCell(cell);
2834  cell->SetUsed(kTRUE);
2835 
2836  }
2837 
2838  } // loop hits last layer
2839  } // loop hits first layer
2840 
2841  } // end seed
2842 
2843 }
2844 
2845 
2846 
2847 
2848 //___________________________________________________________________________
2849 Bool_t AliMFTTrackFinder::LinFit(Int_t nDet, Double_t *xcl,
2850  Double_t *ycl, Double_t *yerr,
2851  Double_t &a, Double_t &ae, Double_t &b, Double_t &be,
2852  Int_t skip) {
2853 
2854 
2855  // y=a*x+b
2856 
2857  const Int_t nMaxh = 100;
2858  Double_t xCl[nMaxh], yCl[nMaxh], yErr[nMaxh];
2859  Int_t idet = 0;
2860  for (Int_t i = 0; i < nDet; i++) {
2861  if (i == skip) continue;
2862  xCl[idet] = xcl[i];
2863  yCl[idet] = ycl[i];
2864  yErr[idet] = yerr[i];
2865  idet++;
2866  }
2867 
2868  Double_t S1, SXY, SX, SY, SXX, SsXY, SsXX, SsYY, Xm, Ym, s, delta, difx;
2869 
2870  S1 = SXY = SX = SY = SXX = 0.0;
2871  SsXX = SsYY = SsXY = Xm = Ym = 0.;
2872  difx = 0.;
2873  for (Int_t i = 0; i < idet; i++) {
2874  S1 += 1.0/(yErr[i]*yErr[i]);
2875  SXY += xCl[i]*yCl[i]/(yErr[i]*yErr[i]);
2876  SX += xCl[i]/(yErr[i]*yErr[i]);
2877  SY += yCl[i]/(yErr[i]*yErr[i]);
2878  SXX += xCl[i]*xCl[i]/(yErr[i]*yErr[i]);
2879  if (i > 0) difx += TMath::Abs(xCl[i]-xCl[i-1]);
2880  Xm += xCl[i];
2881  Ym += yCl[i];
2882  SsXX += xCl[i]*xCl[i];
2883  SsYY += yCl[i]*yCl[i];
2884  SsXY += xCl[i]*yCl[i];
2885  }
2886  delta = SXX*S1 - SX*SX;
2887  if (delta == 0.) {
2888  return kFALSE;
2889  }
2890  a = (SXY*S1 - SX*SY)/delta;
2891  b = (SY*SXX - SX*SXY)/delta;
2892 
2893  Ym /= (Double_t)idet;
2894  Xm /= (Double_t)idet;
2895  SsYY -= (Double_t)idet*(Ym*Ym);
2896  SsXX -= (Double_t)idet*(Xm*Xm);
2897  SsXY -= (Double_t)idet*(Ym*Xm);
2898  Double_t eps = 1.E-24;
2899  if ((idet > 2) && (TMath::Abs(difx) > eps) && ((SsYY-(SsXY*SsXY)/SsXX) > 0.)) {
2900  s = TMath::Sqrt((SsYY-(SsXY*SsXY)/SsXX)/(idet-2));
2901  be = s*TMath::Sqrt(1./(Double_t)idet+(Xm*Xm)/SsXX);
2902  ae = s/TMath::Sqrt(SsXX);
2903  } else {
2904  be = 0.;
2905  ae = 0.;
2906  }
2907 
2908  return kTRUE;
2909 
2910 }
2911 
void CreateCells(Bool_t calcVertex=kFALSE)
Double_t fPlanesZ[fNDetMax]
void CreateCellsOld(Bool_t calcVertex=kFALSE)
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void AddRightNeighbour(Int_t cellgid)
Definition: AliMFTCACell.h:48
const Int_t GetLength()
Definition: AliMFTCACell.h:59
void SetLength(Int_t l1, Int_t l2)
Definition: AliMFTCARoad.h:29
const Int_t GetStatus()
Definition: AliMFTCACell.h:37
void SetInRoad(Int_t i)
Definition: AliMFTCAHit.h:36
void RunBackwardR(AliMFTCARoad *road, Int_t &trackGID)
virtual void Clear(Option_t *)
const Int_t GetNNbL()
Definition: AliMFTCACell.h:49
void DrawTracks(Double_t *pTot, Double_t *Theta)
const Int_t GetNtracks()
void SetChiSqX(Double_t chisq)
Definition: AliMFTCATrack.h:44
const Int_t GetGID() const
Definition: AliMFTCATrack.h:25
Double_t fZVertCalc
! Calculated vertez z
void SetIsMerged()
Definition: AliMFTCACell.h:62
AliMFTCACell * GetCellByGID(Int_t gid)
const Int_t GetDetElemID()
Definition: AliMFTCAHit.h:30
TVectorD vec
Definition: AnalyzeLaser.C:8
AliMFTCAHit * AddHit()
void SetStatus(Int_t s)
Definition: AliMFTCACell.h:20
const Int_t GetLayer1()
Definition: AliMFTCARoad.h:37
const Bool_t IsUsed()
Definition: AliMFTCACell.h:46
void SetID(Int_t id)
Definition: AliMFTCARoad.h:21
Int_t fCellGID
! Cell global identifier
void SetStartLayer(Int_t sl)
Definition: AliMFTCATrack.h:31
Float_t fXCut
Cut in x difference; RuleSelect.
TFile f("CalibObjects.root")
void AddCell(AliMFTCACell *cell)
const Int_t GetNDet()
Double_t fErrX
! Error in X
Int_t * GetLayers()
Definition: AliMFTCACell.h:29
Double_t fErrY
! Error in Y
TClonesArray * fRoads
! Array of roads
void SetDebug(Int_t debug)
AliMFTCALayer * fLayers[fNDetMax]
! Array of layers
AliMFTCALayer * GetLayer(Int_t nl)
const Int_t GetMFTClsId()
Definition: AliMFTCAHit.h:41
Double_t fZVertRange[2]
! Limits of vertex z accepted range
const Int_t GetNcellsInLayer(Int_t nl)
Definition: AliMFTCARoad.h:40
virtual void Clear(const Option_t *)
const Int_t GetNcells()
Definition: AliMFTCALayer.h:30
Double_t * GetHit1()
Definition: AliMFTCACell.h:25
static const Int_t fNMaxPlanes
TH1F * hAngleCells
! Angle between adjacent cells
Double_t GetCellInterceptY(AliMFTCACell *cell, Double_t z)
Double_t fZVertDet
! Vertex z given by ext detector
const Int_t GetNcells() const
Definition: AliMFTCATrack.h:26
Double_t GetY() const
Definition: AliMFTCluster.h:42
void SetUsed()
Definition: AliMFTCAHit.h:33
void PrintCell(Option_t *opt)
Definition: AliMFTCACell.h:65
void ReadParam(Char_t *parfile="param.txt")
void SetDetElemID(Int_t id1, Int_t id2)
Definition: AliMFTCACell.h:24
void SetPhi(Double_t phi)
Definition: AliMFTCATrack.h:38
Double_t * GetHitp1()
Definition: AliMFTCACell.h:27
AliTPCfastTrack * track
const Bool_t IsMerged()
Definition: AliMFTCACell.h:63
const Bool_t HasNbL()
Definition: AliMFTCACell.h:33
TH2F * hDXY[fNDetMax]
! Histogram with X,Y distance between end of cells
void SetMCindex(Int_t index)
Definition: AliMFTCATrack.h:49
void IncrStatus()
Definition: AliMFTCACell.h:36
void PrintTrack(Int_t id)
void sum()
void SetMCflag(UChar_t mcf)
Definition: AliMFTCATrack.h:32
AliMFTCATrack * GetTrack(Int_t nt)
const Int_t GetLayer()
Definition: AliMFTCAHit.h:28
TH1F * hTrackType
! Histogram track Type: 0 = good track ; 1 = hits from other track ; 2 = hits from noisy pixels ...
Int_t fNRoads
! number of built roads
void CreateCellsR(AliMFTCARoad *road)
const Int_t GetTrackGID()
Definition: AliMFTCAHit.h:27
const Int_t GetNbLgid(Int_t id)
Definition: AliMFTCACell.h:51
const Bool_t IsUsed()
Definition: AliMFTCAHit.h:34
Double_t * GetHit2()
Definition: AliMFTCACell.h:26
const Int_t GetNNbR()
Definition: AliMFTCACell.h:50
Float_t fACutV[fNDetMax]
Cut in angle difference: for cell vertex compatibility.
Float_t fYCut
Cut in y difference; RuleSelect.
const Bool_t ReadGeom()
void LoadClusters(TClonesArray *clusterArrayFront[AliMFTConstants::fNMaxPlanes], TClonesArray *clusterArrayBack[AliMFTConstants::fNMaxPlanes])
void SetMFTClsId(Int_t id)
Definition: AliMFTCAHit.h:40
Double_t GetCellInterceptX(AliMFTCACell *cell, Double_t z)
void Reset()
Definition: AliMFTCACell.h:53
void AddLeftNeighbour(Int_t cellgid)
Definition: AliMFTCACell.h:47
AliMFTCATrack * AddTrack(Int_t gid)
void SetGID(Int_t gid, Int_t trackid1, Int_t trackid2)
Definition: AliMFTCACell.h:18
AliMFTCARoad * AddRoad()
Bool_t LinFit(Int_t nDet, Double_t *xcl, Double_t *ycl, Double_t *yerr, Double_t &a, Double_t &ae, Double_t &b, Double_t &be, Int_t skip=-1)
Int_t fNlayers
Number of detection planes.
const Int_t GetNhits()
Definition: AliMFTCALayer.h:29
void SetLayers(Int_t iL1, Int_t iL2)
Definition: AliMFTCACell.h:22
void SetVertZ(Double_t z)
Definition: AliMFTCATrack.h:36
Float_t fACutN[fNDetMax]
Cut in angle difference: for neighbor cells compatibility.
Double_t GetX() const
Definition: AliMFTCluster.h:41
AliMFTCACell * GetCell(Int_t ic) const
Definition: AliMFTCATrack.h:28
Int_t fNtracks
! Number of tracks
Bool_t RuleSelectCell(AliMFTCACell *cell)
Float_t fMaxSegAngle
Max cut of the Theta angle of segments [deg].
Double_t GetCellAngleDif(AliMFTCACell *cell1, AliMFTCACell *cell2)
TH1F * hThetaCells
! Theta of the cells
const Double_t * GetPos()
Definition: AliMFTCAHit.h:19
void SetNHitSta(Int_t n)
Definition: AliMFTCARoad.h:27
AliMFTCACell * GetCellInLayer(Int_t nl, Int_t nc)
Definition: AliMFTCARoad.h:41
const Int_t GetID()
Definition: AliMFTCARoad.h:22
const Int_t GetNcells()
const Bool_t IsGood()
Definition: AliMFTCARoad.h:36
TClonesArray * fTracks
! Array of tracks
const Int_t GetLayer2()
Definition: AliMFTCARoad.h:38
void Init(Char_t *parfile)
const Int_t GetStartLayer() const
Definition: AliMFTCATrack.h:30
Double_t AddCellToChiSq(AliMFTCACell *cell)
TH1F * hNGoodCell
! Histogram showing numbers of good cells in the track
const Int_t GetTrackGID(Int_t p)
Definition: AliMFTCACell.h:39
Double_t fPlaneDetEff[fNDetMax]
TH1F * hDAv[fNDetMax]
! Histogram with angle with respect to the vertex
void SetChiSqY(Double_t chisq)
Definition: AliMFTCATrack.h:45
void SetVertY(Double_t y)
Definition: AliMFTCATrack.h:35
AliMFTCACell * GetCell(Int_t nc)
Definition: AliMFTCALayer.h:32
const Bool_t HasNbR()
Definition: AliMFTCACell.h:34
void AddHit(AliMFTCAHit *hit)
Int_t fMaxCellStatus
! Maximum value of a cell status after RunForward
const Int_t GetCellGID(Int_t ic) const
Definition: AliMFTCATrack.h:27
void SetPos(Double_t x, Double_t y, Double_t z)
Definition: AliMFTCAHit.h:17
const Int_t GetLastCellGID() const
Definition: AliMFTCATrack.h:29
Double_t * GetHitp2()
Definition: AliMFTCACell.h:28
Double_t GetZ() const
Definition: AliMFTCluster.h:43
void RunForwardR(AliMFTCARoad *road, Int_t &trackGID)
const Int_t GetID()
Definition: AliMFTCAHit.h:29
void SetHits(Double_t *h1, Double_t *h2, Double_t z1, Double_t z2)
void SetUsed(Bool_t isu)
Definition: AliMFTCACell.h:45
AliMFTCACell * AddCell()
void AddCell(AliMFTCACell *cell)
AliMFTCAHit * GetHitInLayer(Int_t nl, Int_t nh)
Definition: AliMFTCARoad.h:26
void SetGood()
Definition: AliMFTCARoad.h:35
void SetGID(Int_t gid)
Definition: AliMFTCATrack.h:24
static const Int_t kNDisks
Number of Disk.
const Int_t GetGID()
Definition: AliMFTCACell.h:38
const Int_t IsFace()
Definition: AliMFTCAHit.h:31
const Int_t GetID()
Definition: AliMFTCALayer.h:28
void SetTrackGID(Int_t gid, Int_t la, Int_t id, Int_t detid)
Definition: AliMFTCAHit.h:20
const Int_t GetNhitsInLayer(Int_t nl)
Definition: AliMFTCARoad.h:25
Bool_t RuleSelect2LayersGap(Int_t iL1, Int_t iL2, Double_t *hit1, Double_t *hit2)
Int_t GetMCLabel(Int_t track) const
Definition: AliMFTCluster.h:63
ClassImp(AliMFTTrackFinder) AliMFTTrackFinder
AliMFTCARoad * GetRoad(Int_t nr)
TH1F * hDA[fNDetMax]
! Histogram with angle between cells
void SetMFTClsId(Int_t id1, Int_t id2)
Definition: AliMFTCACell.h:23
TVector3 * GetSeg()
Definition: AliMFTCACell.h:32
AliMFTCACell * GetCellByGID(Int_t gid)
void UpdateStatus()
Definition: AliMFTCACell.h:35
void SetTheta(Double_t the)
Definition: AliMFTCATrack.h:37
void SetVertX(Double_t x)
Definition: AliMFTCATrack.h:34
Bool_t RuleSelect(AliMFTCACell *cellL, AliMFTCACell *cellR)
AliMFTCAHit * GetHit(Int_t nh)
Definition: AliMFTCALayer.h:34