Line data Source code
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 :
41 12 : ClassImp(AliMFTTrackFinder)
42 :
43 : //___________________________________________________________________________
44 0 : AliMFTTrackFinder::AliMFTTrackFinder() :
45 0 : TObject(),
46 0 : fXCut(0.0025),
47 0 : fYCut(0.0025),
48 0 : fMaxSegAngle(180.-15.),
49 0 : fCellGID(0),
50 0 : fMaxCellStatus(0),
51 0 : fNlayers(fNDetMax),
52 0 : fNtracks(0),
53 0 : fCalcVertex(kFALSE),
54 0 : fZVertCalc(0.),
55 0 : fZVertDet(0.),
56 0 : fZVertRange(),
57 0 : fMinTrackLength(fNDetMax),
58 0 : fPlaneDetEff(),
59 0 : fAddNoise(kTRUE),
60 0 : fPixelNoise(1.E-5),
61 0 : fAddQED(kFALSE),
62 0 : fMBRate(0.),
63 0 : fCMOSIntTime(0.),
64 0 : fThick(0.004),
65 0 : fReadGeom(kTRUE),
66 0 : fUseTF(kFALSE),
67 0 : fDebug(0),
68 0 : fCPUTime(0.),
69 0 : fRealTime(0.),
70 0 : fNDifTracks(0),
71 0 : fPlanesZ(),
72 0 : fZGap(1.0),
73 0 : fNRoads(0),
74 0 : fErrX(0.),
75 0 : fErrY(0.)
76 0 : {
77 :
78 : #ifdef OLDGEOM
79 : fZGap = 0.2;
80 : #else
81 0 : fZGap = 1.0;
82 : #endif
83 :
84 0 : fZVertRange[0] = fZVertRange[1] = 0.;
85 :
86 0 : fHistList = new TList();
87 :
88 0 : fTracks = new TClonesArray("AliMFTCATrack", 1000);
89 :
90 0 : fRoads = new TClonesArray("AliMFTCARoad", 1000);
91 :
92 0 : for (Int_t iL = 0; iL < fNlayers; iL++) {
93 :
94 0 : fACutV[iL] = 0.20;
95 0 : fACutN[iL] = 0.10;
96 :
97 0 : fLayers[iL] = new AliMFTCALayer();
98 0 : fLayers[iL]->SetID(iL);
99 0 : fPlaneDetEff[iL] = 1.00;
100 0 : fPlanesZ[iL] = 0.;
101 :
102 0 : 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 0 : 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 0 : fHistList->Add(hDA[iL]);
105 0 : fHistList->Add(hDAv[iL]);
106 0 : 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 0 : fHistList->Add(hDXY[iL]);
109 :
110 : }
111 0 : hNGoodCell = new TH1F("hNGoodCell","Number of good cell in the track",5,-0.5,4.5);
112 0 : fHistList->Add(hNGoodCell);
113 0 : hTrackType = new TH1F("hTrackType","Track type",4,-0.5,3.5);
114 0 : fHistList->Add(hTrackType);
115 0 : hAngleCells = new TH1F("hAngleCell1gap2gaps","#Delta #theta between two segments one 1 gap and one 2 layers gap; #Delta #theta (deg)",200,0.,45.);
116 0 : fHistList->Add(hAngleCells);
117 0 : hThetaCells = new TH1F("hThetaCells","#theta of the cells; #theta (deg)",200,0.,15.);
118 0 : fHistList->Add(hThetaCells);
119 :
120 : // limited range for possible z vertex
121 0 : hVertZ = new TH1F("hVertZ","hVertZ",400,-10.,+10.);
122 0 : hVertZa = new TH1F("hVertZa","hVertZa",400,-10.,+10.);
123 :
124 0 : fHistList->Add(hVertZ);
125 0 : fHistList->Add(hVertZa);
126 :
127 : // temporary
128 0 : hHitDifXY = new TH2F("hHitDifXY","hHitDifXY",100,-0.04,+0.04,100,-0.04,+0.04);
129 0 : fHistList->Add(hHitDifXY);
130 0 : hAngDifAll = new TH1F("hAngDifAll","hAngDifAll",1000,0.,+0.5);
131 0 : fHistList->Add(hAngDifAll);
132 0 : hAngDifDup = new TH1F("hAngDifDup","hAngDifDup",1000,0.,+0.5);
133 0 : fHistList->Add(hAngDifDup);
134 0 : hIntDifXYAll = new TH2F("hIntDifXYAll","hIntDifXYAll",100,-0.01,+0.01,100,-0.01,+0.01);
135 0 : 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 0 : hIntDifXYDup = new TH2F("hIntDifXYDup","hIntDifXYDup",100,-0.1,+0.1,100,-0.1,+0.1);
140 : #endif
141 0 : fHistList->Add(hIntDifXYDup);
142 :
143 0 : }
144 :
145 : //___________________________________________________________________________
146 : AliMFTTrackFinder::~AliMFTTrackFinder()
147 0 : {
148 :
149 0 : }
150 : //___________________________________________________________________________
151 : void AliMFTTrackFinder::Init(Char_t *parfile)
152 : {
153 :
154 0 : AliInfo(Form("Parameter file set to %s ",parfile));
155 :
156 0 : ReadParam(parfile);
157 :
158 0 : ReadGeom();
159 :
160 0 : SetDebug(1);
161 :
162 :
163 :
164 0 : }
165 : //___________________________________________________________________________
166 : void AliMFTTrackFinder::LoadClusters( TClonesArray *clusterArrayFront[AliMFTConstants::fNMaxPlanes], TClonesArray *clusterArrayBack[AliMFTConstants::fNMaxPlanes] ){
167 0 : AliCodeTimerAuto("",0);
168 : AliMFTCALayer *caLayer = NULL;
169 : AliMFTCAHit *caHit = NULL;
170 :
171 0 : for (int iPlane = 0; iPlane<AliMFTConstants::kNDisks; iPlane++) {
172 0 : Int_t nClusterFront = clusterArrayFront[iPlane]->GetEntriesFast();
173 0 : Int_t nClusterBack = clusterArrayBack[iPlane]->GetEntriesFast();
174 0 : 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 0 : caLayer = GetLayer(iPlane*2);
178 0 : for (Int_t iCluster=0; iCluster<nClusterFront; iCluster++) {
179 0 : caHit = caLayer->AddHit();
180 0 : AliMFTCluster * cluster = (AliMFTCluster *)clusterArrayFront[iPlane]->At(iCluster);
181 0 : caHit->SetPos(cluster->GetX(), cluster->GetY(), cluster->GetZ());
182 0 : caHit->SetTrackGID(cluster->GetMCLabel(0),iPlane*2,caLayer->GetNhits()-1,0);
183 0 : caHit->SetMFTClsId(iCluster);
184 : }
185 : // Treating BACK face of the plane
186 0 : caLayer = GetLayer(iPlane*2+1);
187 0 : for (Int_t iCluster=0; iCluster<nClusterBack; iCluster++) {
188 0 : caHit = caLayer->AddHit();
189 0 : AliMFTCluster * cluster = (AliMFTCluster *)clusterArrayBack[iPlane]->At(iCluster);
190 0 : caHit->SetPos(cluster->GetX(), cluster->GetY(), cluster->GetZ());
191 0 : caHit->SetTrackGID(cluster->GetMCLabel(0),iPlane*2+1,caLayer->GetNhits()-1,0);
192 0 : caHit->SetMFTClsId(iCluster);
193 : }
194 :
195 : }
196 0 : }
197 :
198 : //___________________________________________________________________________
199 : void AliMFTTrackFinder::ReadParam(Char_t *parfile)
200 : {
201 :
202 0 : AliInfo(Form("Reading Parameter File %s ",parfile));
203 :
204 0 : std::ifstream in;
205 0 : in.open(parfile,std::ios::in);
206 :
207 0 : std::string line;
208 :
209 0 : getline(in,line);
210 0 : in >> fUseTF;
211 0 : printf("Use the TrackFinder: %d \n",fUseTF);
212 :
213 0 : getline(in,line); getline(in,line);
214 0 : in >> fNlayers;
215 0 : printf("Number of detecting planes: %d \n",fNlayers);
216 :
217 0 : getline(in,line); getline(in,line);
218 0 : in >> fThick;
219 0 : printf("Layer thickness in X0: %5.3f \n",fThick);
220 :
221 0 : getline(in,line); getline(in,line);
222 0 : in >> fPixelNoise;
223 0 : printf("Pixel noise: %4.2e \n",fPixelNoise);
224 0 : in >> fAddNoise;
225 0 : printf("Add noise: %d \n",fAddNoise);
226 :
227 0 : getline(in,line); getline(in,line);
228 0 : in >> fMinTrackLength;
229 0 : printf("fMinTrackLength: %d \n",fMinTrackLength);
230 :
231 0 : getline(in,line); getline(in,line);
232 0 : in >> fXCut;
233 0 : printf("fXCut: %6.4f cm\n",fXCut);
234 0 : in >> fYCut;
235 0 : printf("fYCut: %6.4f cm\n",fYCut);
236 :
237 0 : getline(in,line); getline(in,line);
238 0 : in >> fMaxSegAngle;
239 0 : printf("fMaxSegAngle: %4.1f deg\n",fMaxSegAngle);
240 0 : fMaxSegAngle = 180. - fMaxSegAngle;
241 :
242 0 : getline(in,line); getline(in,line);
243 0 : for (Int_t i = 0; i < fNlayers; i++) {
244 0 : in >> fACutV[i];
245 0 : printf("fACutV[%d]: %4.2f \n",i,fACutV[i]);
246 : }
247 :
248 0 : getline(in,line); getline(in,line);
249 0 : for (Int_t i = 0; i < fNlayers; i++) {
250 0 : in >> fACutN[i];
251 0 : printf("fACutN[%d]: %4.2f \n",i,fACutN[i]);
252 : }
253 :
254 0 : getline(in,line); getline(in,line);
255 0 : for (Int_t i = 0; i < fNlayers; i++) {
256 0 : in >> fPlaneDetEff[i];
257 0 : printf("PlaneDetEff[%d]: %4.2f \n",i,fPlaneDetEff[i]);
258 : }
259 :
260 0 : getline(in,line); getline(in,line);
261 0 : in >> fCalcVertex;
262 0 : printf("Calculate vertex: %d \n",fCalcVertex);
263 :
264 0 : getline(in,line); getline(in,line);
265 0 : in >> fAddQED;
266 0 : printf("Add QED hits: %d \n",fAddQED);
267 0 : in >> fMBRate;
268 0 : printf("Hadronic MB Rate: %5.1f kHz\n",fMBRate);
269 0 : in >> fCMOSIntTime;
270 0 : printf("CMOS integration time: %5.1f microsec \n",fCMOSIntTime);
271 :
272 0 : getline(in,line); getline(in,line);
273 0 : in >> fReadGeom;
274 0 : printf("Read geometry: %d \n",fReadGeom);
275 0 : in >> fGeomName;
276 0 : if (fReadGeom) printf("... from file: %s \n",fGeomName.Data());
277 :
278 0 : in.close();
279 0 : printf("==================================\n");
280 :
281 0 : }
282 :
283 : //___________________________________________________________________________
284 : void AliMFTTrackFinder::Clear(Option_t *)
285 : {
286 :
287 0 : if (fTracks) fTracks->Clear("C");
288 :
289 0 : for (Int_t iL = 0; iL < fNlayers; iL++) {
290 0 : fLayers[iL]->Clear("");
291 : }
292 :
293 0 : fCellGID = 0;
294 0 : fMaxCellStatus = 0;
295 0 : fNtracks = 0;
296 :
297 0 : if (fRoads) fRoads->Clear("C");
298 :
299 0 : fNRoads = 0;
300 :
301 0 : }
302 :
303 : //___________________________________________________________________________
304 : void AliMFTTrackFinder::ClearCells()
305 : {
306 :
307 0 : for (Int_t iL = 0; iL < fNlayers; iL++) {
308 0 : fLayers[iL]->ClearCells();
309 : }
310 :
311 0 : fCellGID = 0;
312 0 : fMaxCellStatus = 0;
313 :
314 0 : }
315 :
316 : //___________________________________________________________________________
317 : void AliMFTTrackFinder::ResetCells()
318 : {
319 :
320 : AliMFTCALayer *layer;
321 : AliMFTCACell *cell;
322 :
323 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
324 0 : layer = GetLayer(iL);
325 0 : for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
326 0 : cell = layer->GetCell(iC);
327 0 : cell->Reset();
328 : }
329 : }
330 :
331 0 : }
332 :
333 : //___________________________________________________________________________
334 : Double_t AliMFTTrackFinder::GetCellAngleDif(AliMFTCACell *cell1, AliMFTCACell *cell2) {
335 :
336 0 : Double_t *c1h1 = cell1->GetHit1();
337 0 : Double_t *c1h2 = cell1->GetHit2();
338 0 : Double_t *c2h1 = cell2->GetHit1();
339 0 : Double_t *c2h2 = cell2->GetHit2();
340 :
341 0 : TVector3 seg1 = TVector3(c1h2[0]-c1h1[0],c1h2[1]-c1h1[1],c1h2[2]-c1h1[2]);
342 0 : TVector3 seg2 = TVector3(c2h2[0]-c2h1[0],c2h2[1]-c2h1[1],c2h2[2]-c2h1[2]);
343 :
344 0 : return (seg1.Angle(seg2))*TMath::RadToDeg();
345 :
346 0 : }
347 :
348 : //___________________________________________________________________________
349 : Double_t AliMFTTrackFinder::GetCellInterceptX(AliMFTCACell *cell, Double_t z) {
350 :
351 0 : Double_t *ch1 = cell->GetHit1();
352 0 : Double_t *ch2 = cell->GetHit2();
353 :
354 0 : return ch1[0]+(ch2[0]-ch1[0])/(ch2[2]-ch1[2])*(z-ch1[2]);
355 :
356 : }
357 :
358 : //___________________________________________________________________________
359 : Double_t AliMFTTrackFinder::GetCellInterceptY(AliMFTCACell *cell, Double_t z) {
360 :
361 0 : Double_t *ch1 = cell->GetHit1();
362 0 : Double_t *ch2 = cell->GetHit2();
363 :
364 0 : return ch1[1]+(ch2[1]-ch1[1])/(ch2[2]-ch1[2])*(z-ch1[2]);
365 :
366 : }
367 :
368 : //___________________________________________________________________________
369 : void AliMFTTrackFinder::AnalyzeCells() {
370 :
371 0 : 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 0 : Double_t hit1[3], hit2[3];
391 : Bool_t multCell1;
392 0 : TList *multCell = new TList();
393 :
394 : const Int_t nMaxh = 100;
395 0 : Double_t xTr[nMaxh], yTr[nMaxh], zTr[nMaxh];
396 0 : Double_t ax, axe, bx, bxe, ay, aye, by, bye;
397 0 : Double_t xTrErrDet = 0.0025/TMath::Sqrt(12.);
398 0 : 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 0 : Double_t xTrErr[nMaxh], yTrErr[nMaxh];
402 0 : for (Int_t i = 0; i < nMaxh; i++) {
403 0 : xTrErr[i] = TMath::Sqrt(xTrErrDet*xTrErrDet+xTrErrMS*xTrErrMS);
404 0 : 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 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
470 :
471 0 : layer = GetLayer(iL);
472 :
473 0 : nCells = layer->GetNcells();
474 :
475 0 : for (Int_t iC1 = 0; iC1 < nCells; iC1++) {
476 :
477 0 : cell1 = layer->GetCell(iC1);
478 0 : if (cell1->GetStatus() < 1 || cell1->IsMerged()) continue;
479 :
480 0 : zc1h1 = cell1->GetHitp1()[2];
481 0 : zc1h2 = cell1->GetHitp2()[2];
482 0 : trkid1h1 = cell1->GetTrackGID(0);
483 0 : trkid1h2 = cell1->GetTrackGID(1);
484 :
485 : multCell1 = kFALSE;
486 0 : for (Int_t iC2 = (iC1+1); iC2 < nCells; iC2++) {
487 :
488 0 : cell2 = layer->GetCell(iC2);
489 0 : if (cell2->GetStatus() < 1 || cell2->IsMerged()) continue;
490 :
491 0 : if (cell2->GetLength() != cell1->GetLength()) continue;
492 :
493 0 : zc2h1 = cell2->GetHitp1()[2];
494 0 : zc2h2 = cell2->GetHitp2()[2];
495 0 : trkid2h1 = cell2->GetTrackGID(0);
496 0 : trkid2h2 = cell2->GetTrackGID(1);
497 :
498 0 : if ((TMath::Abs(zc1h1-zc2h1) > 0.5*fZGap) ||
499 0 : (TMath::Abs(zc1h2-zc2h2) > 0.5*fZGap)) {
500 :
501 0 : xc1h1 = cell1->GetHit1()[0];
502 0 : xc1h2 = cell1->GetHit2()[0];
503 0 : xc2h1 = cell2->GetHit1()[0];
504 0 : xc2h2 = cell2->GetHit2()[0];
505 0 : yc1h1 = cell1->GetHit1()[1];
506 0 : yc1h2 = cell1->GetHit2()[1];
507 0 : yc2h1 = cell2->GetHit1()[1];
508 0 : 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 0 : cellDifX1 = xc1h1-xc2h1;
518 0 : cellDifX2 = xc1h2-xc2h2;
519 0 : cellDifY1 = yc1h1-yc2h1;
520 0 : cellDifY2 = yc1h2-yc2h2;
521 :
522 0 : if (trkid1h1 == trkid2h1 && trkid1h2 == trkid2h2) {
523 0 : cellAngDif = GetCellAngleDif(cell1,cell2);
524 0 : hAngDifDup->Fill(cellAngDif);
525 0 : hIntDifXYDup->Fill(cellDifX1,cellDifY1);
526 0 : hIntDifXYDup->Fill(cellDifX2,cellDifY2);
527 : //printf("CellIntDif %f %f %f %f \n",cellDifX1,cellDifX2,cellDifY1,cellDifY2);
528 0 : }
529 :
530 : //if ((cellAngDif < cellAngDifCut) &&
531 : // (TMath::Abs(cellIntDifX) < fXCut) &&
532 : // (TMath::Abs(cellIntDifY) < fYCut)) {
533 0 : if (TMath::Abs(cellDifX1) < cellIntDifCut &&
534 0 : TMath::Abs(cellDifX2) < cellIntDifCut &&
535 0 : TMath::Abs(cellDifY1) < cellIntDifCut &&
536 0 : TMath::Abs(cellDifY2) < cellIntDifCut) {
537 0 : if (!multCell1) {
538 : //printf("Add mult cell1 %d \n",cell1->GetGID());
539 : //cell1->PrintCell("MC");
540 : multCell1 = kTRUE;
541 0 : multCell->Add(cell1);
542 0 : }
543 : //printf("Add mult cell2 %d \n",cell2->GetGID());
544 : //cell2->PrintCell("MC");
545 0 : multCell->Add(cell2);
546 0 : }
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 0 : 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 0 : for (Int_t imc = 0; imc < multCell->GetSize(); imc++) {
572 0 : cellM = (AliMFTCACell*)multCell->At(imc);
573 0 : xTr[nTr] = cellM->GetHit1()[0];
574 0 : yTr[nTr] = cellM->GetHit1()[1];
575 0 : zTr[nTr] = cellM->GetHit1()[2];
576 0 : zMin = TMath::Min(zMin,zTr[nTr]);
577 0 : zMax = TMath::Max(zMax,zTr[nTr]);
578 0 : nTr++;
579 0 : xTr[nTr] = cellM->GetHit2()[0];
580 0 : yTr[nTr] = cellM->GetHit2()[1];
581 0 : zTr[nTr] = cellM->GetHit2()[2];
582 0 : zMin = TMath::Min(zMin,zTr[nTr]);
583 0 : zMax = TMath::Max(zMax,zTr[nTr]);
584 0 : nTr++;
585 0 : 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 0 : if (LinFit(nTr,zTr,xTr,xTrErr,ax,axe,bx,bxe) &&
596 0 : LinFit(nTr,zTr,yTr,yTrErr,ay,aye,by,bye)) {
597 0 : hit1[0] = ax*zMin+bx;
598 0 : hit1[1] = ay*zMin+by;
599 0 : hit1[2] = zMin;
600 0 : hit2[0] = ax*zMax+bx;
601 0 : hit2[1] = ay*zMax+by;
602 0 : hit2[2] = zMax;
603 : // add a new cell
604 0 : cell = layer->AddCell();
605 0 : cell->SetHits(hit1,hit2,fPlanesZ[cell1->GetLayers()[0]],fPlanesZ[cell1->GetLayers()[1]]);
606 0 : cell->SetLayers(cell1->GetLayers()[0],cell1->GetLayers()[1]);
607 0 : cell->SetStatus(1);
608 0 : cell->SetGID(fCellGID++,trkid1h1,trkid1h2);
609 0 : 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 0 : } else {
621 0 : printf("No line fit possible!\n");
622 : }
623 0 : multCell->Clear();
624 0 : } // end merge multi cells
625 :
626 : } // end cell1 loop
627 :
628 : } // end layer loop
629 :
630 0 : delete multCell;
631 :
632 0 : }
633 :
634 : //___________________________________________________________________________
635 : void AliMFTTrackFinder::CreateGapCells() {
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 0 : Double_t h1[3], h2[3];
645 : Int_t iL2, nH2, trackID1, trackID2;
646 :
647 : Int_t nGapCells = 0;
648 :
649 0 : for (Int_t iL1 = 0; iL1 < (fNlayers-1); iL1++) {
650 :
651 0 : layerL = GetLayer(iL1);
652 :
653 0 : for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) { // Loop over cell in layer iL
654 :
655 : // short cell
656 0 : cellL = layerL->GetCell(iCL);
657 :
658 : // attach at right a long cell
659 0 : if ((iL1 < (fNlayers-3)) && (cellL->GetNNbR() == 0)) {
660 :
661 0 : h1[0] = cellL->GetHit2()[0];
662 0 : h1[1] = cellL->GetHit2()[1];
663 0 : h1[2] = cellL->GetHit2()[2];
664 :
665 0 : iL2 = iL1 + 3;
666 0 : nH2 = GetLayer(iL2)->GetNhits();
667 :
668 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
669 :
670 0 : h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
671 0 : h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
672 0 : h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
673 :
674 0 : TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
675 0 : if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
676 :
677 0 : if (!RuleSelect2LayersGap(iL1+1,iL1+2,h1,h2)) continue;
678 :
679 0 : if (!RuleSelectCell(h1,h2,iL1+1)) continue;
680 :
681 0 : nGapCells++;
682 :
683 0 : AliMFTCACell *cell = GetLayer(iL1+1)->AddCell();
684 0 : cell->SetHits(h1,h2,fPlanesZ[iL1+1],fPlanesZ[iL1+3]);
685 0 : cell->SetLayers(iL1+1,iL1+3);
686 0 : cell->SetStatus(1);
687 0 : trackID1 = cellL->GetTrackGID(1);
688 0 : trackID2 = GetLayer(iL1+3)->GetHit(iH2)->GetTrackGID();
689 0 : cell->SetGID(fCellGID++,trackID1,trackID2);
690 :
691 0 : if (prn) printf("Create gap (L) cell %d in layer %d. \n",cell->GetGID(),iL1+1);
692 :
693 0 : } // end loop hits
694 :
695 0 : } // end attach at right a long cell
696 :
697 : // attach at left a long cell; 1 <> 2
698 0 : if ((iL1 > 1) && (cellL->GetNNbL() == 0)) {
699 :
700 0 : h1[0] = cellL->GetHit1()[0];
701 0 : h1[1] = cellL->GetHit1()[1];
702 0 : h1[2] = cellL->GetHit1()[2];
703 :
704 0 : iL2 = iL1 - 2;
705 0 : nH2 = GetLayer(iL2)->GetNhits();
706 :
707 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
708 :
709 0 : h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
710 0 : h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
711 0 : h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
712 :
713 0 : TVector3 vec(h1[0]-h2[0],h1[1]-h2[1],h1[2]-h2[2]);
714 0 : if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
715 :
716 0 : if (!RuleSelect2LayersGap(iL1-2,iL1-1,h2,h1)) continue;
717 :
718 0 : if (!RuleSelectCell(h2,h1,iL1-2)) continue;
719 :
720 0 : nGapCells++;
721 :
722 0 : AliMFTCACell *cell = GetLayer(iL1-2)->AddCell();
723 0 : cell->SetHits(h2,h1,fPlanesZ[iL1-2],fPlanesZ[iL1]);
724 0 : cell->SetLayers(iL1-2,iL1);
725 0 : cell->SetStatus(1);
726 0 : trackID2 = cellL->GetTrackGID(0);
727 0 : trackID1 = GetLayer(iL1-2)->GetHit(iH2)->GetTrackGID();
728 0 : cell->SetGID(fCellGID++,trackID1,trackID2);
729 :
730 0 : if (prn) printf("Create gap (R) cell %d in layer %d. \n",cell->GetGID(),iL1-2);
731 :
732 0 : } // end loop hits
733 :
734 0 : } // end attach at left a long cell
735 :
736 : } // end loop cells
737 :
738 : } // end loop layers
739 :
740 0 : printf("Found %d gap cells.\n",nGapCells);
741 :
742 0 : }
743 :
744 : //___________________________________________________________________________
745 : void AliMFTTrackFinder::CreateCellsR(AliMFTCARoad *road) {
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 0 : Double_t h1[3], h2[3], h[3], hx, hy, dR;
759 : Int_t mftClsId1, mftClsId2;
760 : Int_t nCombi = 0;
761 :
762 0 : iL1min = road->GetLayer1();
763 0 : iL1max = road->GetLayer2()-1;
764 : //printf("R%d iL1 %d %d \n",ir,iL1min,iL1max);
765 :
766 0 : for (iL1 = iL1min; iL1 <= iL1max; iL1++) {
767 : //printf("iL1 %d \n",iL1);
768 0 : iL2min = iL1+1;
769 0 : nH1 = road->GetNhitsInLayer(iL1);
770 0 : for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
771 : //printf("iH1 %d \n",iH1);
772 0 : hit1 = road->GetHitInLayer(iL1,iH1);
773 0 : iL2max = TMath::Min((iL1+(5-hit1->IsFace())),fNlayers-1);
774 0 : h1[0] = hit1->GetPos()[0];
775 0 : h1[1] = hit1->GetPos()[1];
776 0 : h1[2] = hit1->GetPos()[2];
777 0 : mftClsId1 = hit1->GetMFTClsId();
778 : noCell = kTRUE;
779 : iL2 = iL2min;
780 0 : while (noCell && (iL2 <= iL2max)) {
781 : //printf("iL2 %d \n",iL2);
782 0 : nH2 = road->GetNhitsInLayer(iL2);
783 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
784 0 : hit2 = road->GetHitInLayer(iL2,iH2);
785 0 : h2[0] = hit2->GetPos()[0];
786 0 : h2[1] = hit2->GetPos()[1];
787 0 : h2[2] = hit2->GetPos()[2];
788 0 : mftClsId2 = hit2->GetMFTClsId();
789 0 : nCombi++;
790 0 : if (RuleSelectCell(h1,h2,iL1)) {
791 : noCell = kFALSE;
792 0 : cell = GetLayer(iL1)->AddCell();
793 0 : cell->SetHits(h1,h2,fPlanesZ[iL1],fPlanesZ[iL2]);
794 0 : cell->SetMFTClsId(mftClsId1,mftClsId2);
795 0 : cell->SetLayers(iL1,iL2);
796 0 : cell->SetStatus(1);
797 0 : trackID1 = hit1->GetTrackGID();
798 0 : trackID2 = hit2->GetTrackGID();
799 0 : detElemID1 = hit1->GetDetElemID();
800 0 : detElemID2 = hit2->GetDetElemID();
801 0 : cell->SetGID(fCellGID++,trackID1,trackID2);
802 0 : cell->SetDetElemID(detElemID1,detElemID2);
803 0 : 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 0 : } // end create cell
806 : } // end loop iH2
807 0 : 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 0 : }
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 0 : Double_t h1[3], h2[3], h[3], hx, hy, dR;
857 : Int_t nCombi = 0;
858 :
859 : iL1min = 0;
860 0 : iL1max = fNlayers-2;
861 :
862 0 : for (iL1 = iL1min; iL1 <= iL1max; iL1++) {
863 : //printf("iL1 %d \n",iL1);
864 0 : iL2min = iL1+1;
865 0 : nH1 = GetLayer(iL1)->GetNhits();
866 0 : for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
867 : //printf("iH1 %d \n",iH1);
868 0 : hit1 = GetLayer(iL1)->GetHit(iH1);
869 0 : iL2max = TMath::Min((iL1+(5-hit1->IsFace())),fNlayers-1);
870 0 : h1[0] = hit1->GetPos()[0];
871 0 : h1[1] = hit1->GetPos()[1];
872 0 : h1[2] = hit1->GetPos()[2];
873 : noCell = kTRUE;
874 : iL2 = iL2min;
875 0 : while (noCell && (iL2 <= iL2max)) {
876 : //printf("iL2 %d \n",iL2);
877 0 : nH2 = GetLayer(iL2)->GetNhits();
878 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
879 0 : hit2 = GetLayer(iL2)->GetHit(iH2);
880 0 : h2[0] = hit2->GetPos()[0];
881 0 : h2[1] = hit2->GetPos()[1];
882 0 : h2[2] = hit2->GetPos()[2];
883 0 : nCombi++;
884 0 : if (RuleSelectCell(h1,h2,iL1)) {
885 : noCell = kFALSE;
886 0 : cell = GetLayer(iL1)->AddCell();
887 0 : cell->SetHits(h1,h2,fPlanesZ[iL1],fPlanesZ[iL2]);
888 0 : cell->SetLayers(iL1,iL2);
889 0 : cell->SetStatus(1);
890 0 : trackID1 = hit1->GetTrackGID();
891 0 : trackID2 = hit2->GetTrackGID();
892 0 : detElemID1 = hit1->GetDetElemID();
893 0 : detElemID2 = hit2->GetDetElemID();
894 0 : cell->SetGID(fCellGID++,trackID1,trackID2);
895 0 : 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 0 : } // end create cell
898 : } // end loop iH2
899 0 : 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 0 : }
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 0 : Double_t h1[3], h2[3];
967 : Double_t nCombi = 0.;
968 :
969 : Int_t iL1, iL2;
970 :
971 : iL1 = 0;
972 0 : while (iL1 < (fNlayers-1)) {
973 :
974 0 : iL2 = iL1 + 1;
975 :
976 0 : nH1 = GetLayer(iL1)->GetNhits();
977 0 : if (prn) printf("---> 1st Layer L %d with %d hits.\n",iL1,nH1);
978 :
979 0 : nH2 = GetLayer(iL2)->GetNhits();
980 0 : if (prn) printf("---> 2nd Layer R %d with %d hits.\n",iL2,nH2);
981 :
982 0 : for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
983 :
984 0 : h1[0] = GetLayer(iL1)->GetHit(iH1)->GetPos()[0];
985 0 : h1[1] = GetLayer(iL1)->GetHit(iH1)->GetPos()[1];
986 0 : h1[2] = GetLayer(iL1)->GetHit(iH1)->GetPos()[2];
987 :
988 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
989 :
990 0 : h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
991 0 : h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
992 0 : h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
993 :
994 0 : TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
995 0 : if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
996 :
997 0 : nCombi += 1.;
998 :
999 0 : if (!cv && !RuleSelectCell(h1,h2,iL1)) continue;
1000 :
1001 0 : AliMFTCACell *cell = GetLayer(iL1)->AddCell();
1002 0 : cell->SetHits(h1,h2,fPlanesZ[iL1],fPlanesZ[iL2]);
1003 0 : cell->SetLayers(iL1,iL2);
1004 0 : cell->SetStatus(1);
1005 0 : trackID1 = GetLayer(iL1)->GetHit(iH1)->GetTrackGID();
1006 0 : trackID2 = GetLayer(iL2)->GetHit(iH2)->GetTrackGID();
1007 0 : detElemID1 = GetLayer(iL1)->GetHit(iH1)->GetDetElemID();
1008 0 : detElemID2 = GetLayer(iL2)->GetHit(iH2)->GetDetElemID();
1009 0 : cell->SetGID(fCellGID++,trackID1,trackID2);
1010 0 : 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 0 : } // end loop hit in layer 2
1017 :
1018 : } // end loop hit in layer 1
1019 :
1020 0 : if (prn) printf("Create cell %d in layer %d-%d. \n",GetLayer(iL1)->GetNcells(),iL1,GetLayer(iL1)->GetID());
1021 :
1022 0 : if (cv) break;
1023 :
1024 : iL1++;
1025 :
1026 : } // end loop layer 1
1027 :
1028 0 : if (cv) CalculateVertex();
1029 :
1030 : if (kTRUE || prn) {
1031 0 : printf("From %.0f combinations: \n",nCombi);
1032 : Long_t nTotCell =0;
1033 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
1034 0 : printf("Layer %d nr of cells: %d \n",iL,GetLayer(iL)->GetNcells());
1035 0 : nTotCell += GetLayer(iL)->GetNcells();
1036 : }
1037 0 : printf("Tot cells: %ld \n",nTotCell);
1038 : }
1039 :
1040 0 : }
1041 :
1042 : //___________________________________________________________________________
1043 : void AliMFTTrackFinder::CalculateVertex() {
1044 :
1045 0 : hVertZ->Reset();
1046 0 : hVertZa->Reset();
1047 :
1048 0 : 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 0 : for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
1056 :
1057 0 : cell = layer->GetCell(iC);
1058 :
1059 0 : h1 = cell->GetHit1();
1060 0 : h2 = cell->GetHit2();
1061 :
1062 0 : x0 = h1[0];
1063 0 : y0 = h1[1];
1064 0 : z0 = h1[2];
1065 0 : a = (h2[0]-h1[0])/(h2[2]-h1[2]);
1066 0 : b = (h2[1]-h1[1])/(h2[2]-h1[2]);
1067 : c = 1.;
1068 :
1069 0 : n1 = (c*x0-a*z0)/c;
1070 0 : n2 = a/c;
1071 0 : n3 = (c*y0-b*z0)/c;
1072 0 : n4 = b/c;
1073 :
1074 0 : zmin = -(n1*n2+n3*n4)/(n2*n2+n4*n4);
1075 :
1076 0 : hVertZ->Fill(zmin);
1077 0 : 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 0 : Int_t binMax = hVertZ->GetNbinsX();
1085 :
1086 0 : hVertZ->Fit("pol1","QW0");
1087 0 : TF1 *f = hVertZ->GetFunction("pol1");
1088 :
1089 0 : for (Int_t i = binMin; i <= binMax; i++) {
1090 0 : hVertZ->SetBinContent(i,TMath::Max(0.,hVertZ->GetBinContent(i)-f->Eval(hVertZ->GetBinCenter(i))));
1091 : }
1092 :
1093 0 : for (Int_t i = binMin; i <= binMax; i++) {
1094 0 : bin1 = TMath::Max(binMin,i-binW);
1095 0 : bin2 = TMath::Min(binMax,i+binW);
1096 0 : sum = hVertZ->Integral(bin1,bin2);
1097 0 : if (sum > maxsum) {
1098 : maxsum = sum;
1099 0 : zvert = hVertZ->GetBinCenter(i);
1100 0 : }
1101 : //printf("%d %d %d %f %f %f \n",bin1,i,bin2,sum,maxsum,zvert);
1102 : }
1103 0 : fZVertCalc = zvert;
1104 0 : printf("Fit vertex z = %f cm \n",fZVertCalc);
1105 :
1106 : // range = zvert +/- 3 cm
1107 0 : fZVertRange[0] = fZVertCalc - 3.;
1108 0 : fZVertRange[1] = fZVertCalc + 3.;
1109 :
1110 0 : }
1111 :
1112 : //___________________________________________________________________________
1113 : void AliMFTTrackFinder::RunForwardR(AliMFTCARoad *road, Int_t& trackGID) {
1114 :
1115 : Bool_t prn = kFALSE;
1116 :
1117 0 : 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 0 : Int_t cellLayers[2];
1128 :
1129 0 : fMaxCellStatus = 0;
1130 :
1131 : iter = 0;
1132 :
1133 : stch = kTRUE;
1134 0 : while (stch) {
1135 :
1136 : stch = kFALSE;
1137 0 : iter++;
1138 :
1139 0 : for (iL = 0; iL < (fNlayers-2); iL++) { // loop over layers
1140 :
1141 0 : for (Int_t iCL = 0; iCL < road->GetNcellsInLayer(iL); iCL++) {
1142 :
1143 0 : cellL = road->GetCellInLayer(iL,iCL);
1144 0 : if (cellL->GetStatus() == 0) continue;
1145 :
1146 0 : for(Int_t i = 0; i < 2; i++) cellLayers[i] = cellL->GetLayers()[i];
1147 :
1148 0 : iR = iL+(cellLayers[1]-cellLayers[0]);
1149 0 : if (iR >= (fNlayers-1))
1150 : continue;
1151 :
1152 0 : for (Int_t iCR = 0; iCR < road->GetNcellsInLayer(iR); iCR++) {
1153 :
1154 0 : cellR = road->GetCellInLayer(iR,iCR);
1155 0 : if (cellR->GetStatus() == 0) continue;
1156 :
1157 0 : nCombiTot += 1;
1158 :
1159 0 : if ((cellL->GetStatus() == cellR->GetStatus()) &&
1160 0 : RuleSelect(cellL,cellR)) {
1161 :
1162 0 : if (prn){
1163 0 : AliInfo(Form("Matching cells: L(%d) cellGID(%d) %d R(%d) cellGID(%d) %d \n",iL,iCL,cellL->GetGID(),iR,iCR,cellR->GetGID()));
1164 0 : }
1165 :
1166 0 : nCombiMatch += 1;
1167 :
1168 0 : if (iter == 1) {
1169 0 : cellL->AddRightNeighbour(cellR->GetGID());
1170 0 : cellR->AddLeftNeighbour(cellL->GetGID());
1171 0 : }
1172 :
1173 0 : cellR->IncrStatus();
1174 :
1175 : stch = kTRUE;
1176 :
1177 0 : } // END : matching cells
1178 :
1179 : } // END : loop over cells in layer iR
1180 :
1181 0 : } // END : loop over cells in layer iL
1182 :
1183 : } // END : loop over layer iL
1184 :
1185 0 : UpdateCellStatusR();
1186 :
1187 0 : if (kFALSE || prn) {
1188 0 : AliInfo(Form("Iteration: %5d ----------------- \n",iter));
1189 0 : for (iL = 0; iL < (fNlayers-1); iL++) {
1190 0 : for (Int_t iCL = 0; iCL < road->GetNcellsInLayer(iL); iCL++) {
1191 0 : cellL = road->GetCellInLayer(iL,iCL);
1192 0 : if (cellL->HasNbL() || cellL->HasNbR()) {
1193 0 : printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1194 0 : }
1195 : }
1196 : }
1197 : }
1198 :
1199 : } // end status change
1200 :
1201 0 : if (kFALSE || prn) {
1202 0 : printf("End iteration: ----------------- \n");
1203 0 : for (iL = 0; iL < (fNlayers-1); iL++) {
1204 0 : for (Int_t iCL = 0; iCL < road->GetNcellsInLayer(iL); iCL++) {
1205 0 : cellL = road->GetCellInLayer(iL,iCL);
1206 0 : if (cellL->HasNbL() || cellL->HasNbR()) {
1207 0 : printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1208 0 : }
1209 : }
1210 : }
1211 : }
1212 :
1213 0 : if (kFALSE || prn) {
1214 :
1215 0 : printf("RunForward after %d iterations, nr of Total combinations %.0f , nr of matched combinations %.0f\n",iter,nCombiTot, nCombiMatch);
1216 :
1217 0 : printf("After RunForward max cell status = %d \n",fMaxCellStatus);
1218 :
1219 0 : }
1220 :
1221 0 : RunBackwardR(road,trackGID);
1222 :
1223 0 : }
1224 :
1225 : //___________________________________________________________________________
1226 : void AliMFTTrackFinder::RunForward() {
1227 :
1228 : Bool_t prn = kFALSE;
1229 :
1230 0 : 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 0 : while (stch) {
1244 :
1245 : nCombiIter = 0;
1246 : nCombiIterMatch = 0;
1247 :
1248 : stch = kFALSE;
1249 0 : iter++;
1250 :
1251 0 : for (iL = 0; iL < (fNlayers-2); iL++) { // loop over layers
1252 :
1253 0 : layerL = GetLayer(iL);
1254 :
1255 0 : if (prn) printf("L %d cells %d \n",iL,layerL->GetNcells());
1256 :
1257 0 : for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) { // Loop over cell in layer iL
1258 :
1259 0 : cellL = layerL->GetCell(iCL);
1260 0 : if (cellL->GetStatus() == 0) continue;
1261 :
1262 0 : Int_t cellLayers[2];
1263 0 : for(Int_t i = 0; i < 2; i++) cellLayers[i] = cellL->GetLayers()[i];
1264 :
1265 0 : iR = iL+(cellLayers[1] - cellLayers[0]);
1266 0 : if (iR < (fNlayers-1))
1267 0 : layerR = GetLayer(iR);
1268 : else
1269 0 : continue;
1270 :
1271 : // vertex selection here ?
1272 : //if (!RuleSelectCell(cellL)) continue;
1273 :
1274 0 : for (Int_t iCR = 0; iCR < layerR->GetNcells(); iCR++) { // Loop over cell in layer iL
1275 :
1276 0 : cellR = layerR->GetCell(iCR);
1277 0 : if (cellR->GetStatus() == 0) continue;
1278 :
1279 : // vertex selection here ?
1280 : //if (!RuleSelectCell(cellR)) continue;
1281 :
1282 0 : nCombiTot += 1;
1283 0 : nCombiIter += 1;
1284 :
1285 0 : if ((cellL->GetStatus() == cellR->GetStatus()) && RuleSelect(cellL,cellR)) { // If Cells are matching in angles and have equal status values
1286 :
1287 0 : if (prn){
1288 0 : printf("Matching cells: L cellGID %d R cellGID %d \n",cellL->GetGID(),cellR->GetGID());
1289 0 : printf("Layer L %d cell %d - Layer R %d cell %d \n",iL,iCL,iR,iCR);
1290 0 : }
1291 0 : nCombiMatch += 1;
1292 0 : nCombiIterMatch += 1;
1293 :
1294 0 : if (iter == 1) {
1295 0 : cellL->AddRightNeighbour(cellR->GetGID());
1296 0 : cellR->AddLeftNeighbour(cellL->GetGID());
1297 0 : }
1298 :
1299 0 : cellR->IncrStatus();
1300 :
1301 : stch = kTRUE;
1302 :
1303 0 : } // END : matching cells
1304 :
1305 : } // END : loop over cell in layer iL-1
1306 :
1307 0 : } // END : loop over cell layer iL
1308 : } // END : loop over layers
1309 :
1310 0 : if (prn) {
1311 0 : printf("Iteration: %5d nr of combinations %.0f match %.0f \n",iter,nCombiIter,nCombiIterMatch);
1312 0 : }
1313 :
1314 0 : UpdateCellStatus();
1315 :
1316 0 : if (prn) {
1317 0 : printf("Iteration: %5d ----------------- \n",iter);
1318 0 : for (iL = 0; iL < (fNlayers-1); iL++) {
1319 0 : layerL = GetLayer(iL);
1320 0 : for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) {
1321 0 : cellL = layerL->GetCell(iCL);
1322 0 : if (cellL->HasNbL() || cellL->HasNbR()) {
1323 0 : printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1324 0 : }
1325 : }
1326 : }
1327 : }
1328 :
1329 : } // end status change
1330 :
1331 0 : if (prn) {
1332 0 : printf("End iteration: ----------------- \n");
1333 0 : for (iL = 0; iL < (fNlayers-1); iL++) {
1334 0 : layerL = GetLayer(iL);
1335 0 : for (Int_t iCL = 0; iCL < layerL->GetNcells(); iCL++) {
1336 0 : cellL = layerL->GetCell(iCL);
1337 0 : if (cellL->HasNbL() || cellL->HasNbR()) {
1338 0 : printf("L%1d C%03d S%1d GID%03d NNb %d %d \n",iL,iCL,cellL->GetStatus(),cellL->GetGID(),cellL->GetNNbL(),cellL->GetNNbR());
1339 0 : }
1340 : }
1341 : }
1342 : }
1343 :
1344 : if (kTRUE || prn) {
1345 :
1346 0 : printf("RunForward after %d iterations, nr of Total combinations %.0f , nr of matched combinations %.0f\n",iter,nCombiTot, nCombiMatch);
1347 :
1348 0 : printf("After RunForward max cell status = %d \n",fMaxCellStatus);
1349 :
1350 : }
1351 :
1352 0 : }
1353 :
1354 : //___________________________________________________________________________
1355 : void AliMFTTrackFinder::RunBackwardR(AliMFTCARoad *road, Int_t& trackGID) {
1356 :
1357 : Bool_t prn = kFALSE;
1358 :
1359 0 : if (prn) printf("Run backward road %d ==================================== \n",road->GetID());
1360 :
1361 0 : 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;
1370 : AliMFTCATrack *track;
1371 :
1372 0 : Bool_t addCellToTrack, hitSta[5];
1373 :
1374 : Int_t minStartLayer = 6;
1375 : Int_t maxStartLayer = 8;
1376 :
1377 0 : for (Int_t startLayer = maxStartLayer; startLayer >= minStartLayer; startLayer--) {
1378 :
1379 0 : if (prn) printf("Start layer %d \n",startLayer);
1380 :
1381 0 : for (Int_t iCR = 0; iCR < road->GetNcellsInLayer(startLayer); iCR++) {
1382 :
1383 0 : cellR = road->GetCellInLayer(startLayer,iCR);
1384 :
1385 0 : if (cellR->GetStatus() == 0) continue;
1386 0 : if (cellR->IsUsed()) continue;
1387 0 : if (cellR->GetStatus() < (fMinTrackLength-1)) continue;
1388 :
1389 0 : if (prn) printf("Create new track %d \n",trackGID);
1390 :
1391 0 : track = AddTrack(trackGID++);
1392 0 : track->SetStartLayer(startLayer);
1393 0 : track->AddCell(cellR);
1394 0 : cellR->SetUsed(kTRUE);
1395 :
1396 : // add cells to new track
1397 : addCellToTrack = kTRUE;
1398 0 : while (addCellToTrack) {
1399 :
1400 0 : 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 0 : for (Int_t iNNbL = 0; iNNbL < cellR->GetNNbL(); iNNbL++) {
1414 :
1415 0 : 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 0 : if (cellL->GetStatus() == 0) continue;
1422 0 : if (cellL->IsUsed()) continue;
1423 0 : if (cellL->GetStatus() != (cellR->GetStatus()-1)) continue;
1424 :
1425 : // ... smallest deviation
1426 0 : cellAngDif = GetCellAngleDif(cellL,cellR);
1427 0 : if (cellAngDifPrev < 0.) {
1428 : cellAngDifPrev = cellAngDif;
1429 0 : } else {
1430 0 : if (cellAngDif < cellAngDifPrev) {
1431 : iCellAngDifMin = iNNbL;
1432 0 : }
1433 : }
1434 :
1435 : // ... smallest chisquare
1436 0 : if (track->AddCellToChiSq(cellL) < chisqNbLprev) {
1437 : iNbLchiSqMin = iNNbL;
1438 0 : }
1439 :
1440 0 : chisqNbLprev = track->AddCellToChiSq(cellL);
1441 :
1442 0 : } // 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 0 : if (cellR->GetNNbL() > 0) {
1453 0 : cellL = road->GetCellByGID(cellR->GetNbLgid(iNbLchiSqMin));
1454 : addCellToTrack = kTRUE;
1455 0 : addCellIdToTrack = cellR->GetNbLgid(iNbLchiSqMin);
1456 0 : cellL = road->GetCellByGID(addCellIdToTrack);
1457 0 : track->AddCell(cellL);
1458 0 : cellL->SetUsed(kTRUE);
1459 0 : if (prn) {
1460 0 : 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 0 : }
1462 : }
1463 :
1464 : } // END : addCellToTrack
1465 :
1466 : // check again track length
1467 0 : for (Int_t j = 0; j < 5; j++) hitSta[j] = kFALSE;
1468 0 : for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1469 0 : cell = GetCellByGID(track->GetCellGID(iCell));
1470 0 : hitSta[cell->GetLayers()[0]/2] = kTRUE;
1471 0 : hitSta[cell->GetLayers()[1]/2] = kTRUE;
1472 : }
1473 : nHitSta = 0;
1474 0 : for (Int_t i = 0; i < 5; i++) {
1475 0 : if (hitSta[i]) nHitSta++;
1476 : }
1477 0 : if (nHitSta < fMinTrackLength) {
1478 0 : for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1479 0 : cell = track->GetCell(iCell);
1480 0 : cell->SetUsed(kFALSE);
1481 : }
1482 0 : RemoveLastTrack();
1483 0 : trackGID--;
1484 0 : } else {
1485 0 : road->SetGood();
1486 : }
1487 :
1488 : } // END : startLayer cells loop
1489 :
1490 : } // END : startLayer loop
1491 :
1492 0 : }
1493 :
1494 : //___________________________________________________________________________
1495 : void AliMFTTrackFinder::RunBackward() {
1496 :
1497 : Bool_t prn = kFALSE;
1498 :
1499 0 : if (prn) printf("Run backward ==================================== \n");
1500 :
1501 0 : 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;
1510 : AliMFTCATrack *track;
1511 :
1512 : Bool_t addCellToTrack;
1513 :
1514 0 : Int_t minStartLayer = fMinTrackLength-2;
1515 0 : Int_t maxStartLayer = (fNlayers-1)-1;
1516 :
1517 0 : for (Int_t startLayer = maxStartLayer; startLayer >= minStartLayer; startLayer--) {
1518 :
1519 0 : if (prn) printf("Start layer %d \n",startLayer);
1520 :
1521 0 : layerR = GetLayer(startLayer);
1522 :
1523 0 : for (Int_t iCR = 0; iCR < layerR->GetNcells(); iCR++) {
1524 0 : cellR = layerR->GetCell(iCR);
1525 0 : if (cellR->GetStatus() == 0) continue;
1526 0 : if (cellR->IsUsed()) continue;
1527 0 : if (cellR->GetStatus() >= (fMinTrackLength-1)) {
1528 0 : if (prn) printf("Create new track %d \n",trackGID);
1529 0 : track = AddTrack(trackGID++);
1530 0 : track->SetStartLayer(startLayer);
1531 : //track->AddCellGID(cellR->GetGID());
1532 0 : track->AddCell(cellR);
1533 0 : cellR->SetUsed(kTRUE);
1534 0 : if (prn) {
1535 0 : 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 0 : }
1537 : // add cells to new track
1538 : addCellToTrack = kTRUE;
1539 0 : while (addCellToTrack) {
1540 0 : cellR = GetCellByGID(track->GetLastCellGID());
1541 : addCellToTrack = kFALSE;
1542 : iNbLchiSqMin = 0;
1543 : chisqNbLprev = 0.;
1544 : cellAngDifPrev = -1.;
1545 : iCellAngDifMin = 0;
1546 0 : for (Int_t iNNbL = 0; iNNbL < cellR->GetNNbL(); iNNbL++) {
1547 0 : cellL = GetCellByGID(cellR->GetNbLgid(iNNbL));
1548 0 : if (cellL->GetStatus() == 0) continue;
1549 0 : if (cellL->IsUsed()) continue;
1550 0 : if (cellL->GetStatus() == (cellR->GetStatus()-1)) {
1551 0 : cellAngDif = GetCellAngleDif(cellL,cellR);
1552 0 : if (cellAngDifPrev < 0.) {
1553 : cellAngDifPrev = cellAngDif;
1554 0 : } else {
1555 0 : if (cellAngDif < cellAngDifPrev) {
1556 : iCellAngDifMin = iNNbL;
1557 0 : }
1558 : }
1559 0 : if (track->AddCellToChiSq(cellL) < chisqNbLprev) {
1560 : iNbLchiSqMin = iNNbL;
1561 0 : }
1562 : }
1563 0 : chisqNbLprev = track->AddCellToChiSq(cellL);
1564 0 : } // 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 0 : if (cellR->GetNNbL() > 0) {
1571 0 : cellL = GetCellByGID(cellR->GetNbLgid(iNbLchiSqMin));
1572 : addCellToTrack = kTRUE;
1573 0 : addCellIdToTrack = cellR->GetNbLgid(iNbLchiSqMin);
1574 : //track->AddCellGID(addCellIdToTrack);
1575 0 : cellL = GetCellByGID(addCellIdToTrack);
1576 0 : track->AddCell(cellL);
1577 0 : cellL->SetUsed(kTRUE);
1578 0 : if (prn) {
1579 0 : 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 0 : }
1581 : }
1582 : } // END : addCellToTrack
1583 :
1584 : // check again track length
1585 0 : if (track->GetNcells() < (fMinTrackLength-1)) {
1586 0 : for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1587 0 : cell = track->GetCell(iCell);
1588 0 : cell->SetUsed(kFALSE);
1589 : }
1590 0 : RemoveLastTrack();
1591 : trackGID--;
1592 0 : }
1593 :
1594 : } // END : create new track
1595 :
1596 : } // END : startLayer cells loop
1597 :
1598 : } // END : startLayer loop
1599 :
1600 0 : }
1601 :
1602 : //___________________________________________________________________________
1603 : void AliMFTTrackFinder::FilterTracks() {
1604 0 : AliCodeTimerAuto("",0);
1605 :
1606 : Bool_t prn = kFALSE;
1607 :
1608 0 : AliInfo(Form("Filtering %d tracks",GetNtracks()));
1609 :
1610 : Int_t nTrkC = 0, nTrkG = 0, nTrkF = 0, nTrkN = 0;
1611 :
1612 : AliMFTCATrack *track;
1613 : AliMFTCACell *cell, *celln;
1614 : Int_t ndof, nptr, cellGID, cellGIDn, cellGIDprev = -1, nTrkSplitEnd = 0;
1615 : const Int_t nMaxh = 100;
1616 0 : Double_t xTr[nMaxh], yTr[nMaxh], zTr[nMaxh];
1617 0 : 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 0 : Int_t nrHitTrackID[nTrackMax], idHitTrackID[nTrackMax], nDiffTracks;
1622 : Int_t idMaxHitsTrack, nMaxHits, nNoisyPix;
1623 0 : for (Int_t i = 0; i < nTrackMax; i++) {
1624 0 : nrHitTrackID[i] = 0;
1625 0 : idHitTrackID[i] = -2;
1626 : }
1627 0 : Int_t nrAliMFTCATrackID[nTrackMax], idAliMFTCATrackID[nTrackMax], nDiffAliMFTCATracks = 0;
1628 0 : for (Int_t i = 0; i < nTrackMax; i++) {
1629 0 : nrAliMFTCATrackID[i] = 0;
1630 0 : idAliMFTCATrackID[i] = -2;
1631 : }
1632 0 : Double_t xTrErrDet = 0.0028/TMath::Sqrt(12.);
1633 0 : 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 0 : Double_t xTrErr[nMaxh], yTrErr[nMaxh];
1637 0 : for (Int_t i = 0; i < nMaxh; i++) {
1638 0 : xTrErr[i] = TMath::Sqrt(xTrErrDet*xTrErrDet+xTrErrMS*xTrErrMS);
1639 0 : yTrErr[i] = TMath::Sqrt(yTrErrDet*yTrErrDet+yTrErrMS*yTrErrMS);
1640 : }
1641 0 : fErrX = xTrErr[0];
1642 0 : fErrY = yTrErr[0];
1643 :
1644 : Int_t nTotalHits = 0, nCleanTotalHits = 0;
1645 :
1646 0 : for (Int_t iTrack = 0; iTrack < GetNtracks(); iTrack++) {
1647 0 : track = GetTrack(iTrack);
1648 : nDiffTracks = 0;
1649 : nptr = 0;
1650 0 : for (Int_t iCell = 0; iCell < track->GetNcells(); iCell++) {
1651 0 : cellGID = track->GetCellGID(iCell);
1652 0 : cell = GetCellByGID(cellGID);
1653 : //printf("Track %3d Cell %5d \n",iTrack,cellGID);
1654 : // tracks split from the first (highest status) cell
1655 0 : if (iCell == 0) {
1656 0 : if (cellGIDprev >= 0) {
1657 0 : if (cellGID == cellGIDprev) {
1658 : // this is a split track
1659 0 : nTrkSplitEnd++;
1660 : splitTrack = kTRUE;
1661 0 : } else {
1662 : splitTrack = kFALSE;
1663 : }
1664 : }
1665 : cellGIDprev = cellGID;
1666 0 : }
1667 0 : if (prn) {
1668 0 : printf("Cell %4d from track %d ",cellGID,track->GetGID());
1669 0 : printf("(%4d %4d %d) \n",cell->GetTrackGID(0),cell->GetTrackGID(1),cell->GetGID());
1670 0 : printf("2: %f %f %f \n",cell->GetHit2()[0],cell->GetHit2()[1],cell->GetHit2()[2]);
1671 0 : printf("1: %f %f %f \n",cell->GetHit1()[0],cell->GetHit1()[1],cell->GetHit1()[2]);
1672 : }
1673 : // extract hit x,y,z
1674 0 : if (nptr == 0) {
1675 0 : xTr[nptr] = cell->GetHit2()[0];
1676 0 : yTr[nptr] = cell->GetHit2()[1];
1677 0 : zTr[nptr] = cell->GetHit2()[2];
1678 0 : nptr++;
1679 0 : xTr[nptr] = cell->GetHit1()[0];
1680 0 : yTr[nptr] = cell->GetHit1()[1];
1681 0 : zTr[nptr] = cell->GetHit1()[2];
1682 0 : nptr++;
1683 0 : } else {
1684 0 : xTr[nptr] = cell->GetHit1()[0];
1685 0 : yTr[nptr] = cell->GetHit1()[1];
1686 0 : zTr[nptr] = cell->GetHit1()[2];
1687 0 : nptr++;
1688 : }
1689 : // count the generated tracks which contribute to this reconstructed track
1690 0 : for (Int_t ihc = 0; ihc < 2; ihc++) {
1691 : recTrack = kTRUE;
1692 0 : for (Int_t idt = 0; idt < nDiffTracks; idt++) {
1693 0 : if (idHitTrackID[idt] == cell->GetTrackGID(ihc)) {
1694 0 : nrHitTrackID[idt]++;
1695 : recTrack = kFALSE;
1696 0 : break;
1697 : }
1698 : }
1699 0 : if (recTrack) {
1700 0 : idHitTrackID[nDiffTracks] = cell->GetTrackGID(ihc);
1701 0 : nrHitTrackID[nDiffTracks] = 1;
1702 0 : nDiffTracks++;
1703 0 : }
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 0 : if (nDiffTracks == 1 && idHitTrackID[0] >= 0) {
1723 : cleanTrack = kTRUE;
1724 : idMaxHitsTrack = idHitTrackID[0];
1725 0 : } else {
1726 : nNoisyPix = 0;
1727 : nMaxHits = 0;
1728 0 : for (Int_t idt = 0; idt < nDiffTracks; idt++) {
1729 0 : if (idHitTrackID[idt] == -1) {
1730 0 : nNoisyPix = nrHitTrackID[idt];
1731 0 : } else if (nMaxHits < nrHitTrackID[idt]) {
1732 : nMaxHits = nrHitTrackID[idt];
1733 : idMaxHitsTrack = idHitTrackID[idt];
1734 0 : }
1735 : }
1736 0 : if (GetNDet() == 5) {
1737 : // allow one fake hit
1738 0 : if (nMaxHits >= (2*track->GetNcells())-1) {
1739 : goodTrack = kTRUE;
1740 0 : } else {
1741 0 : if (nNoisyPix > 0) noiseTrack = kTRUE;
1742 : else fakeTrack = kTRUE;
1743 : }
1744 : }
1745 0 : if (GetNDet() == 5*2) {
1746 : // allow two fake hits
1747 0 : if (nMaxHits >= (2*track->GetNcells())-2) {
1748 : goodTrack = kTRUE;
1749 0 : } else {
1750 0 : if (nNoisyPix > 0) noiseTrack = kTRUE;
1751 : else fakeTrack = kTRUE;
1752 : }
1753 : }
1754 0 : if (GetNDet() == 6) {
1755 : // allow one fake hit
1756 0 : if (nMaxHits >= (2*track->GetNcells())-1) {
1757 : goodTrack = kTRUE;
1758 0 : } else {
1759 0 : if (nNoisyPix > 0) noiseTrack = kTRUE;
1760 : else fakeTrack = kTRUE;
1761 : }
1762 : }
1763 : } // end assert track quality
1764 0 : nTotalHits += nptr;
1765 0 : if (cleanTrack) {
1766 : // count the duplicated clean tracks
1767 : recTrack = kTRUE;
1768 0 : for (Int_t idt = 0; idt < nDiffAliMFTCATracks; idt++) {
1769 0 : if (idAliMFTCATrackID[idt] == idMaxHitsTrack) {
1770 0 : nrAliMFTCATrackID[idt]++;
1771 : recTrack = kFALSE;
1772 0 : break;
1773 : }
1774 : }
1775 0 : if (recTrack) {
1776 0 : idAliMFTCATrackID[nDiffAliMFTCATracks] = idMaxHitsTrack;
1777 0 : nrAliMFTCATrackID[nDiffAliMFTCATracks] = 1;
1778 0 : nDiffAliMFTCATracks++;
1779 0 : }
1780 0 : track->SetMCflag(1);
1781 0 : track->SetMCindex(idMaxHitsTrack);
1782 0 : nTrkC++;
1783 0 : nCleanTotalHits += nptr;
1784 0 : }
1785 0 : if (goodTrack) {
1786 0 : track->SetMCflag(2);
1787 0 : nTrkG++;
1788 0 : }
1789 0 : if (fakeTrack) {
1790 0 : track->SetMCflag(3);
1791 0 : nTrkF++;
1792 0 : }
1793 0 : if (noiseTrack) {
1794 0 : track->SetMCflag(4);
1795 0 : nTrkN++;
1796 0 : }
1797 : // linear regression
1798 : //printf("Fit line with %d points.\n",nptr);
1799 0 : if (LinFit(nptr,zTr,xTr,xTrErr,a,ae,b,be)) {
1800 0 : x0 = b; xS = a;
1801 0 : if (LinFit(nptr,zTr,yTr,yTrErr,a,ae,b,be)) {
1802 0 : y0 = b; yS = a;
1803 : chisqx = 0.;
1804 : chisqy = 0.;
1805 0 : for (Int_t iptr = 0; iptr < nptr; iptr++) {
1806 : //printf("%d %f %f %f \n",iptr,xTr[iptr],yTr[iptr],zTr[iptr]);
1807 0 : chisqx += (xTr[iptr]-(xS*zTr[iptr]+x0))*(xTr[iptr]-(xS*zTr[iptr]+x0))/(xTrErr[iptr]*xTrErr[iptr]);
1808 0 : 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 0 : if (TMath::Abs(xS) > 0.) {
1813 0 : trkPhi = TMath::ATan(yS/xS);
1814 : // put the correct signs
1815 0 : if (xS < 0. && yS > 0.) {
1816 0 : trkPhi = TMath::Pi()+trkPhi;
1817 0 : }
1818 0 : if (xS < 0. && yS < 0.) {
1819 0 : trkPhi = TMath::Pi()+trkPhi;
1820 0 : }
1821 0 : if (xS > 0. && yS < 0.) {
1822 0 : trkPhi = TMath::TwoPi()+trkPhi;
1823 0 : }
1824 0 : if (TMath::Abs(TMath::Sin(trkPhi)) > 0.) {
1825 0 : trkThe = TMath::ATan(yS/TMath::Sin(trkPhi));
1826 : // put the correct signs
1827 0 : trkThe = TMath::Abs(trkThe); // is always smaller than 90deg
1828 0 : trkPhi *= TMath::RadToDeg();
1829 0 : trkThe *= TMath::RadToDeg();
1830 0 : }
1831 : }
1832 : // calculate DCA with the beam axis
1833 0 : zmin = -(x0*xS+y0*yS)/(xS*xS+yS*yS);
1834 0 : track->SetTheta(trkThe);
1835 0 : track->SetPhi(trkPhi);
1836 0 : if (fCalcVertex) {
1837 0 : track->SetVertX(x0+xS*fZVertCalc);
1838 0 : track->SetVertY(y0+yS*fZVertCalc);
1839 0 : track->SetVertZ(fZVertCalc);
1840 0 : } else {
1841 0 : track->SetVertX(x0+xS*fZVertDet);
1842 0 : track->SetVertY(y0+yS*fZVertDet);
1843 0 : track->SetVertZ(fZVertDet);
1844 : }
1845 0 : ndof = nptr-2;
1846 0 : track->SetChiSqX(chisqx/(Double_t)ndof);
1847 0 : track->SetChiSqY(chisqy/(Double_t)ndof);
1848 0 : } // yz fit
1849 : } // xz fit
1850 : //printf("End fit.\n");
1851 : } // end track loop
1852 :
1853 0 : fNDifTracks = nDiffAliMFTCATracks;
1854 :
1855 0 : 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 0 : }
1858 :
1859 : //___________________________________________________________________________
1860 : void AliMFTTrackFinder::DrawTracks(Double_t *pTot, Double_t *Theta) {
1861 :
1862 : Bool_t prn = kFALSE;
1863 :
1864 : AliMFTCATrack *track;
1865 : AliMFTCACell *cell;
1866 : Int_t cellGID, trackGID, nTrackS = 0;
1867 0 : Bool_t single[10000], hitFromNoisyPix, hitFromDiffTrack;
1868 0 : Int_t nGoodCell, firstHitTrackID, idTrack[10000], nGoodTracks = 0;
1869 0 : for (Int_t i = 0; i < 10000; i++) idTrack[i] = -1;
1870 :
1871 0 : 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 0 : for (Int_t iT = 0; iT < GetNtracks(); iT++) {
1882 : // cell content in a track
1883 0 : single[iT] = kTRUE;
1884 0 : track = GetTrack(iT);
1885 0 : cellGID = track->GetCellGID(0);
1886 0 : cell = GetCellByGID(cellGID);
1887 0 : trackGID = cell->GetTrackGID(0);
1888 0 : single[iT] &= (cell->GetTrackGID(1) == trackGID);
1889 : nGoodCell = 0;
1890 0 : for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1891 0 : cellGID = track->GetCellGID(iC);
1892 0 : cell = GetCellByGID(cellGID);
1893 0 : single[iT] &= (cell->GetTrackGID(0) == trackGID);
1894 0 : single[iT] &= (cell->GetTrackGID(1) == trackGID);
1895 0 : 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 0 : if (fDebug > 0) {
1923 0 : hNGoodCell->Fill(nGoodCell);
1924 0 : }
1925 0 : if (single[iT]) {
1926 0 : nTrackS++;
1927 0 : } else {
1928 0 : if (prn || fDebug > 0) {
1929 0 : for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1930 0 : cellGID = track->GetCellGID(iC);
1931 0 : cell = GetCellByGID(cellGID);
1932 0 : if (prn) printf("Track %4d Cell %2d TrackGID %5d %5d \n",iT,iC,cell->GetTrackGID(0),cell->GetTrackGID(1));
1933 : }
1934 0 : }
1935 : }
1936 0 : if (prn) printf("Track %d with %d cells: \n",iT,track->GetNcells());
1937 : hitFromNoisyPix = kFALSE;
1938 : hitFromDiffTrack = kFALSE;
1939 : firstHitTrackID = -1;
1940 0 : for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
1941 0 : cellGID = track->GetCellGID(iC);
1942 0 : cell = GetCellByGID(cellGID);
1943 0 : if (iC == 0) firstHitTrackID = cell->GetTrackGID(0);
1944 0 : if (prn) printf("\t%5d from track: ",cellGID);
1945 0 : if (prn) printf("%5d %5d \n",cell->GetTrackGID(0),cell->GetTrackGID(1));
1946 0 : if (cell->GetTrackGID(0) == -1 || cell->GetTrackGID(1) == -1)
1947 0 : hitFromNoisyPix = kTRUE;
1948 0 : else if (cell->GetTrackGID(0) != firstHitTrackID ||
1949 0 : cell->GetTrackGID(1) != firstHitTrackID)
1950 0 : hitFromDiffTrack = kTRUE;
1951 :
1952 : //cell->PrintCell("FULL");
1953 : }
1954 0 : if (hitFromDiffTrack) {
1955 0 : hTrackType->Fill(1);
1956 : //PrintTrack(iT);
1957 0 : } else if (hitFromNoisyPix) {
1958 0 : hTrackType->Fill(2);
1959 : //PrintTrack(iT);
1960 0 : } else {
1961 0 : if (idTrack[firstHitTrackID] >= 0) {
1962 0 : hTrackType->Fill(3);
1963 0 : printf("Double track %6d from %6d and %6d \n",firstHitTrackID,iT,idTrack[firstHitTrackID]);
1964 0 : PrintTrack(iT);
1965 0 : PrintTrack(idTrack[firstHitTrackID]);
1966 0 : } else {
1967 0 : idTrack[firstHitTrackID] = iT;
1968 0 : hTrackType->Fill(0);
1969 0 : if (track->GetNcells() >= (fMinTrackLength-1)) nGoodTracks++;
1970 : //PrintTrack(iT);
1971 : //printf("%d \n",track->GetNcells());
1972 : }
1973 : }
1974 : } // end loop tracks
1975 :
1976 0 : printf("... %d single. \n",nTrackS);
1977 0 : printf("... %d (%d) good \n",(Int_t)hTrackType->GetBinContent(1),nGoodTracks);
1978 :
1979 0 : }
1980 :
1981 : //___________________________________________________________________________
1982 : AliMFTCACell *AliMFTTrackFinder::GetCellByGID(Int_t gid) {
1983 :
1984 : AliMFTCALayer *layer;
1985 : AliMFTCACell *cell;
1986 :
1987 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
1988 0 : layer = GetLayer(iL);
1989 0 : for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
1990 0 : cell = layer->GetCell(iC);
1991 0 : if (gid == cell->GetGID()) return cell;
1992 : }
1993 : }
1994 :
1995 0 : return 0;
1996 :
1997 0 : }
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 0 : Int_t layer = cellL->GetLayers()[1];
2013 :
2014 0 : TVector3 *segL_ = cellL->GetSeg();
2015 0 : TVector3 *segR_ = cellR->GetSeg();
2016 :
2017 0 : TVector3 segL = TVector3(*segL_);
2018 0 : TVector3 segR = TVector3(*segR_);
2019 :
2020 : const Double_t *hitL[2];
2021 : const Double_t *hitR[2];
2022 :
2023 0 : hitL[0] = cellL->GetHit1();
2024 0 : hitL[1] = cellL->GetHit2();
2025 0 : hitR[0] = cellR->GetHit1();
2026 0 : 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 0 : 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 0 : 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 0 : dx = (dx > 0.) ? (TMath::Sqrt(dx)) : 0.;
2034 0 : dy = (dy > 0.) ? (TMath::Sqrt(dy)) : 0.;
2035 0 : a = (segL.Angle(segR))*TMath::RadToDeg(); // Angle between the segments of each cell
2036 :
2037 0 : if (fDebug > 0) {
2038 0 : hDA[cellL->GetLayers()[1]]->Fill(a);
2039 0 : 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 0 : if ((dx > fXCut) || (dy > fYCut) || (a > fACutN[layer])) return kFALSE;
2073 :
2074 0 : return kTRUE;
2075 :
2076 0 : }
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 0 : if (iL1 >= 0) layer1 = GetLayer(iL1);
2087 0 : else return kFALSE;
2088 :
2089 : AliMFTCALayer *layer2 = 0;
2090 0 : if (iL2 >= 0) layer2 = GetLayer(iL2);
2091 :
2092 0 : for (Int_t iCell = 0; iCell < layer1->GetNcells(); iCell++) {
2093 :
2094 0 : AliMFTCACell * cell = layer1->GetCell(iCell);
2095 0 : 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 0 : 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 0 : 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 0 : if ((TMath::Abs(dx) > fXCut) || (TMath::Abs(dy) > fYCut)) continue;
2106 :
2107 0 : TVector3 *seg1_ = cell->GetSeg();
2108 0 : TVector3 seg1 = TVector3(*seg1_);
2109 :
2110 0 : TVector3 seg2(hit2[0]-hit1[0],hit2[1]-hit1[1],hit2[2]-hit1[2]);
2111 0 : Double_t a = (seg1.Angle(seg2))*TMath::RadToDeg(); // Angle between the segments of each cell
2112 :
2113 : //printf(Form("Angle = %f\n",a));
2114 0 : if (fDebug > 0) hAngleCells->Fill(a);
2115 :
2116 0 : if (prn) {
2117 0 : printf("Compare with cell %d in layer %d \n",iCell,layer1->GetID());
2118 : }
2119 :
2120 0 : if (a < 0.1) findCompatibleCells = kTRUE;
2121 :
2122 0 : }
2123 :
2124 0 : if (iL2 < 0) return (!findCompatibleCells);
2125 :
2126 0 : for (Int_t iCell = 0; iCell < layer2->GetNcells(); iCell++) {
2127 :
2128 0 : AliMFTCACell * cell = layer2->GetCell(iCell);
2129 0 : 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 0 : 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 0 : 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 0 : if ((TMath::Abs(dx) > fXCut) || (TMath::Abs(dy) > fYCut)) continue;
2140 :
2141 0 : TVector3 *seg1_ = cell->GetSeg();
2142 0 : TVector3 seg1 = TVector3(*seg1_);
2143 :
2144 0 : TVector3 seg2(hit2[0]-hit1[0],hit2[1]-hit1[1],hit2[2]-hit1[2]);
2145 0 : Double_t a = (seg1.Angle(seg2))*TMath::RadToDeg(); // Angle between the segments of each cell
2146 :
2147 : //printf(Form("Angle = %f\n",a));
2148 0 : if (fDebug > 0) hAngleCells->Fill(a);
2149 :
2150 0 : if (prn) {
2151 0 : printf("Compare with cell %d in layer %d \n",iCell,layer2->GetID());
2152 : }
2153 :
2154 0 : if (a < 0.1) findCompatibleCells = kTRUE;
2155 :
2156 0 : }
2157 :
2158 0 : return (!findCompatibleCells);
2159 :
2160 0 : }
2161 :
2162 :
2163 :
2164 : //___________________________________________________________________________
2165 : Bool_t AliMFTTrackFinder::RuleSelectCell(AliMFTCACell *cell) { // Look if segment pointing to the vertex
2166 :
2167 0 : Int_t layer = cell->GetLayers()[0];
2168 :
2169 0 : TVector3 *seg_ = cell->GetSeg();
2170 :
2171 0 : TVector3 seg = TVector3(*seg_);
2172 :
2173 : const Double_t *hit;
2174 :
2175 0 : hit = cell->GetHit1();
2176 :
2177 : Double_t vert[3] = { 0., 0., 0. };
2178 0 : if (fCalcVertex) vert[2] = fZVertCalc;
2179 0 : else vert[2] = fZVertDet;
2180 :
2181 : Double_t a;
2182 0 : TVector3 segV;
2183 :
2184 0 : segV = TVector3(hit[0]-vert[0],hit[1]-vert[1],hit[2]-vert[2]);
2185 0 : a = (seg.Angle(segV))*TMath::RadToDeg();
2186 :
2187 0 : hThetaCells->Fill(seg.Theta()*TMath::RadToDeg());
2188 0 : hDAv[layer]->Fill(a);
2189 :
2190 0 : if ( a > fACutV[layer]) return kFALSE;
2191 :
2192 0 : return kTRUE;
2193 :
2194 0 : }
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 0 : TVector3 segV, segVdet;
2203 :
2204 0 : 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 0 : if (fCalcVertex) vert[2] = fZVertCalc;
2208 0 : else vert[2] = fZVertDet;
2209 :
2210 0 : segVdet = TVector3(h1[0]-vert[0],h1[1]-vert[1],h1[2]-vert[2]);
2211 0 : a = (seg.Angle(segVdet))*TMath::RadToDeg();
2212 :
2213 : //hDAv[iL1]->Fill(a);
2214 :
2215 0 : if (acalc) {
2216 0 : acut = f->Eval(180.-segVdet.Theta()*TMath::RadToDeg());
2217 0 : } else {
2218 0 : 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 0 : if ( a > acut) {
2224 0 : AliDebug(2,"Cell NOT Selected");
2225 0 : return kFALSE;
2226 : }
2227 0 : AliDebug(2,"Cell Selected");
2228 :
2229 0 : return kTRUE;
2230 :
2231 0 : }
2232 :
2233 : //___________________________________________________________________________
2234 : AliMFTCATrack* AliMFTTrackFinder::AddTrack(Int_t gid) {
2235 :
2236 0 : new ((*fTracks)[fNtracks++]) AliMFTCATrack();
2237 0 : AliMFTCATrack *track = (AliMFTCATrack*)fTracks->At(fTracks->GetLast());
2238 0 : track->SetGID(gid);
2239 :
2240 0 : return track;
2241 :
2242 0 : }
2243 :
2244 : //___________________________________________________________________________
2245 : AliMFTCATrack* AliMFTTrackFinder::AddTrack(Int_t gid, const AliMFTCATrack& trk) {
2246 :
2247 : // create new track and copy
2248 :
2249 0 : new ((*fTracks)[fNtracks++]) AliMFTCATrack();
2250 0 : AliMFTCATrack *track = (AliMFTCATrack*)fTracks->At(fTracks->GetLast());
2251 0 : track->SetGID(gid);
2252 :
2253 0 : track->SetStartLayer(trk.GetStartLayer());
2254 :
2255 0 : for (Int_t iC = 0; iC < trk.GetNcells(); iC++) {
2256 : //track->AddCellGID(trk.GetCellGID(iC));
2257 0 : track->AddCell(trk.GetCell(iC));
2258 : }
2259 :
2260 0 : return track;
2261 :
2262 0 : }
2263 :
2264 : //___________________________________________________________________________
2265 : void AliMFTTrackFinder::UpdateCellStatusR() {
2266 :
2267 : AliMFTCARoad *road;
2268 : AliMFTCACell *cell;
2269 :
2270 0 : for (Int_t ir = 0; ir < fNRoads; ir++) {
2271 0 : road = GetRoad(ir);
2272 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
2273 0 : for (Int_t iC = 0; iC < road->GetNcellsInLayer(iL); iC++) {
2274 0 : cell = road->GetCellInLayer(iL,iC);
2275 0 : cell->UpdateStatus();
2276 0 : fMaxCellStatus = TMath::Max(fMaxCellStatus,cell->GetStatus());
2277 : }
2278 : }
2279 : }
2280 :
2281 0 : }
2282 :
2283 : //___________________________________________________________________________
2284 : void AliMFTTrackFinder::UpdateCellStatus() {
2285 :
2286 : AliMFTCALayer *layer;
2287 : AliMFTCACell *cell;
2288 :
2289 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
2290 0 : layer = GetLayer(iL);
2291 0 : for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
2292 0 : cell = layer->GetCell(iC);
2293 0 : cell->UpdateStatus();
2294 0 : fMaxCellStatus = TMath::Max(fMaxCellStatus,cell->GetStatus());
2295 : }
2296 : }
2297 :
2298 0 : }
2299 :
2300 : //___________________________________________________________________________
2301 : void AliMFTTrackFinder::PrintTrack(Int_t id) {
2302 :
2303 0 : AliMFTCATrack *track = GetTrack(id);
2304 : Int_t cellGID;
2305 : AliMFTCACell *cell;
2306 0 : printf("Track:\t%6d \n",id);
2307 0 : for (Int_t iC = 0; iC < track->GetNcells(); iC++) {
2308 0 : cellGID = track->GetCellGID(iC);
2309 0 : cell = GetCellByGID(cellGID);
2310 0 : printf("cell\t%6d gid\t%6d \tfrom track: ",iC,cellGID);
2311 0 : printf("\t%6d\t%6d \tin layers: ",cell->GetTrackGID(0),cell->GetTrackGID(1));
2312 0 : printf("\t%1d\t%1d\n",cell->GetLayers()[0],cell->GetLayers()[1]);
2313 : }
2314 :
2315 0 : }
2316 :
2317 : //___________________________________________________________________________
2318 : void AliMFTTrackFinder::PrintAll() {
2319 :
2320 : AliMFTCALayer *layer;
2321 : AliMFTCACell *cell;
2322 :
2323 0 : for (Int_t iL = 0; iL < (fNlayers-1); iL++) {
2324 0 : layer = GetLayer(iL);
2325 0 : printf("LayerID %d \n",layer->GetID());
2326 0 : for (Int_t iC = 0; iC < layer->GetNcells(); iC++) {
2327 0 : cell = layer->GetCell(iC);
2328 0 : cell->PrintCell("FULL");
2329 : }
2330 : }
2331 :
2332 0 : }
2333 :
2334 : //___________________________________________________________________________
2335 : void AliMFTTrackFinder::DrawHisto() {
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 0 : }
2355 :
2356 : //___________________________________________________________________________
2357 : void AliMFTTrackFinder::PrintParam() {
2358 :
2359 0 : printf("==== AliMFTTrackFinder::PrintParam ====\n");
2360 :
2361 0 : printf("Number of layers: %d \n",fNlayers);
2362 0 : printf("Use the TrackFinder: %d \n",fUseTF);
2363 0 : printf("Layer thickness in X0: %5.3f \n",fThick);
2364 0 : printf("Pixel noise: %4.2e \n",fPixelNoise);
2365 0 : printf("Add noise: %d \n",fAddNoise);
2366 0 : printf("fMinTrackLength: %d \n",fMinTrackLength);
2367 0 : printf("fXCut: %6.4f cm\n",fXCut);
2368 0 : printf("fYCut: %6.4f cm\n",fYCut);
2369 0 : printf("fMaxSegAngle: %4.1f deg\n",fMaxSegAngle);
2370 0 : for (Int_t i = 0; i < fNlayers; i++) {
2371 0 : printf("fACutV[%d]: %4.2f \n",i,fACutV[i]);
2372 : }
2373 0 : for (Int_t i = 0; i < fNlayers; i++) {
2374 0 : printf("fACutN[%d]: %4.2f \n",i,fACutN[i]);
2375 : }
2376 0 : for (Int_t i = 0; i < fNlayers; i++) {
2377 0 : printf("PlaneDetEff[%d]: %4.2f \n",i,fPlaneDetEff[i]);
2378 : }
2379 0 : printf("Calculate vertex: %d \n",fCalcVertex);
2380 :
2381 0 : printf("CPU time: %f seconds\n",fCPUTime);
2382 0 : printf("Real time: %f seconds\n",fRealTime);
2383 :
2384 0 : printf("============================\n");
2385 :
2386 0 : }
2387 :
2388 : //___________________________________________________________________________
2389 : AliMFTCARoad* AliMFTTrackFinder::AddRoad() {
2390 :
2391 0 : new ((*fRoads)[fNRoads++]) AliMFTCARoad();
2392 0 : AliMFTCARoad *road = (AliMFTCARoad*)fRoads->At(fRoads->GetLast());
2393 :
2394 0 : return road;
2395 :
2396 0 : }
2397 :
2398 : //___________________________________________________________________________
2399 : void AliMFTTrackFinder::BuildRoads() {
2400 0 : 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 0 : Int_t iPla2Min[4] = { 6, 6, 6, 6 };
2413 0 : Int_t iPla2Max[4] = { 9, 9, 7, 7 };
2414 :
2415 : Int_t nH1, nH2, nH;
2416 0 : Int_t roadLen, nRoads = 0, trackGID = 0;
2417 0 : 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 0 : Bool_t hitSta[5] = { kFALSE, kFALSE, kFALSE, kFALSE, kFALSE };
2424 : Bool_t isUsed, newRoad;
2425 :
2426 0 : for (Int_t iL1 = iPla1Min; iL1 <= iPla1Max; iL1++) {
2427 :
2428 0 : for (Int_t iL2 = iPla2Max[iL1]; iL2 >= iPla2Min[iL1]; iL2--) {
2429 :
2430 : // see AliMFTCARoad::SetLength
2431 0 : roadLen = iL2/2-iL1/2;
2432 :
2433 0 : nH1 = GetLayer(iL1)->GetNhits();
2434 0 : if(prn) AliInfo(Form("nH1 = %d ",nH1));
2435 0 : for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
2436 : //AliInfo(Form("iH1 = %d ",iH1));
2437 :
2438 0 : hit1 = GetLayer(iL1)->GetHit(iH1);
2439 :
2440 0 : 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 0 : if (prn)
2456 0 : printf("Hit1: %d %d %d %d \n",hit1->GetLayer(),hit1->GetID(),hit1->GetDetElemID(),hit1->GetTrackGID());
2457 :
2458 0 : h1[0] = hit1->GetPos()[0];
2459 0 : h1[1] = hit1->GetPos()[1];
2460 0 : h1[2] = hit1->GetPos()[2];
2461 :
2462 0 : nH2 = GetLayer(iL2)->GetNhits();
2463 0 : if(prn) AliInfo(Form("nH2 = %d ",nH2));
2464 :
2465 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
2466 :
2467 0 : hit2 = GetLayer(iL2)->GetHit(iH2);
2468 :
2469 0 : 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 0 : if (prn)
2485 0 : printf("Hit2: %d %d %d %d \n",hit2->GetLayer(),hit2->GetID(),hit2->GetDetElemID(),hit2->GetTrackGID());
2486 :
2487 0 : h2[0] = hit2->GetPos()[0];
2488 0 : h2[1] = hit2->GetPos()[1];
2489 0 : h2[2] = hit2->GetPos()[2];
2490 :
2491 0 : TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
2492 0 : if (vec.Theta() < fMaxSegAngle*TMath::DegToRad()) continue;
2493 :
2494 0 : if (!RuleSelectCell(h1,h2,iL1)) continue;
2495 :
2496 : newRoad = kTRUE;
2497 0 : for (Int_t i = 0; i < 5; i++) hitSta[i] = kFALSE;
2498 :
2499 0 : for (Int_t iL = (iL1+1); iL <= (iL2-1); iL++) {
2500 :
2501 0 : nH = GetLayer(iL)->GetNhits();
2502 :
2503 0 : for (Int_t iH = 0; iH < nH; iH++) {
2504 :
2505 0 : hit = GetLayer(iL)->GetHit(iH);
2506 :
2507 0 : if (hit->IsUsed()) continue;
2508 :
2509 0 : h[0] = hit->GetPos()[0];
2510 0 : h[1] = hit->GetPos()[1];
2511 0 : h[2] = hit->GetPos()[2];
2512 :
2513 0 : hx = h1[0] + (h2[0]-h1[0])*(h[2]-h1[2])/(h2[2]-h1[2]);
2514 0 : hy = h1[1] + (h2[1]-h1[1])*(h[2]-h1[2])/(h2[2]-h1[2]);
2515 :
2516 0 : dR = TMath::Sqrt((hx-h[0])*(hx-h[0])+(hy-h[1])*(hy-h[1]));
2517 0 : 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 0 : if (prn)
2533 0 : printf("Hit: %d %d %d %d \n",hit->GetLayer(),hit->GetID(),hit->GetDetElemID(),hit->GetTrackGID());
2534 :
2535 0 : if (newRoad) {
2536 : //AliInfo("Adding new road");
2537 0 : road = AddRoad();
2538 0 : road->SetID(nRoads++);
2539 0 : road->AddHit(hit1);
2540 0 : road->AddHit(hit2);
2541 :
2542 0 : hit1->SetInRoad(road->GetID());
2543 0 : hit2->SetInRoad(road->GetID());
2544 :
2545 0 : hitSta[iL1/2] = kTRUE;
2546 0 : hitSta[iL2/2] = kTRUE;
2547 :
2548 0 : road->SetLength(iL1,iL2);
2549 :
2550 : newRoad = kFALSE;
2551 :
2552 0 : }
2553 :
2554 0 : road->AddHit(hit);
2555 0 : hit->SetInRoad(road->GetID());
2556 0 : hitSta[iL/2] = kTRUE;
2557 :
2558 0 : } // end loop hits middle plane
2559 :
2560 : } // end loop middle plane
2561 :
2562 : // count the number of hit stations (disks)
2563 0 : if (!newRoad) {
2564 : Int_t nHitSta = 0;
2565 0 : for (Int_t i = 0; i < 5; i++)
2566 0 : if (hitSta[i]) nHitSta++;
2567 0 : road->SetNHitSta(nHitSta);
2568 0 : if (nHitSta >= fMinTrackLength) {
2569 0 : CreateCellsR(road);
2570 0 : RunForwardR(road,trackGID);
2571 0 : if (road->IsGood()) {
2572 0 : for (Int_t j = 0; j < fNlayers; j++) {
2573 0 : for (Int_t k = 0; k < road->GetNhitsInLayer(j); k++) {
2574 0 : hit = road->GetHitInLayer(j,k);
2575 : //printf("%d %d %d \n",hit->GetLayer(),hit->GetID(),GetLayer(hit->GetLayer())->GetNhits());
2576 0 : htmp = GetLayer(hit->GetLayer())->GetHit(hit->GetID());
2577 0 : htmp->SetUsed();
2578 0 : if (prn)
2579 0 : printf("Hit used: %d %d %d %d \n",hit->GetLayer(),hit->GetID(),hit->GetDetElemID(),hit->GetTrackGID());
2580 : }
2581 : }
2582 0 : }
2583 : }
2584 0 : }
2585 :
2586 0 : } // end loop hits last plane
2587 :
2588 0 : } // 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 0 : }
2606 :
2607 : //___________________________________________________________________________
2608 : void AliMFTTrackFinder::FindTracks() {
2609 0 : AliCodeTimerAuto("",0);
2610 :
2611 : Bool_t prn = kFALSE;
2612 :
2613 0 : Double_t h1[3], h2[3], h[3], hx, hy;
2614 0 : Double_t htr1[3], htr2[3];
2615 : const Int_t nMaxh = 100;
2616 0 : Double_t trX[nMaxh], trY[nMaxh], trZ[nMaxh];
2617 0 : 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;
2623 : AliMFTCATrack *track;
2624 :
2625 : Double_t dR, dRmin, dRcut = 0.0100;
2626 :
2627 0 : 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 0 : par[1] = (cut177-cut170)/(177.-170.);
2632 0 : par[0] = cut170-par[1]*170.;
2633 0 : fACutF->SetParameter(0,par[0]);
2634 0 : fACutF->SetParameter(1,par[1]);
2635 :
2636 : Bool_t addHit, selCAgapCell;
2637 0 : Bool_t hitSta[5];
2638 : Bool_t seed = kTRUE;
2639 :
2640 : Int_t step = 0;
2641 :
2642 : iL1 = 0;
2643 :
2644 0 : while (seed) {
2645 :
2646 0 : if (step == 0) {
2647 0 : iL2 = fNlayers-1;
2648 0 : } else {
2649 0 : iL2--;
2650 : }
2651 :
2652 0 : step++;
2653 :
2654 0 : if (iL2 < iL1 + (fMinTrackLength-1)) {
2655 0 : iL1++;
2656 0 : if (iL1 > fNlayers - (fMinTrackLength-1)) break;
2657 : step = 0;
2658 0 : continue;
2659 : }
2660 :
2661 0 : if (prn) printf("iL1,iL2 %d %d \n",iL1,iL2);
2662 : //continue;
2663 :
2664 0 : nH1 = GetLayer(iL1)->GetNhits();
2665 0 : nH2 = GetLayer(iL2)->GetNhits();
2666 0 : if (prn) printf("nH1 %d nH2 %d\n",nH1,nH2);
2667 0 : for (Int_t iH1 = 0; iH1 < nH1; iH1++) {
2668 :
2669 0 : if (GetLayer(iL1)->GetHit(iH1)->IsUsed()) continue;
2670 :
2671 :
2672 0 : h1[0] = GetLayer(iL1)->GetHit(iH1)->GetPos()[0];
2673 0 : h1[1] = GetLayer(iL1)->GetHit(iH1)->GetPos()[1];
2674 0 : h1[2] = GetLayer(iL1)->GetHit(iH1)->GetPos()[2];
2675 :
2676 0 : if (prn) printf("H1 %d (%.1f,%.1f,%.1f)\n",iH1,h1[0],h1[1],h1[2]);
2677 :
2678 0 : for (Int_t iH2 = 0; iH2 < nH2; iH2++) {
2679 :
2680 0 : if (GetLayer(iL2)->GetHit(iH2)->IsUsed()) continue;
2681 :
2682 :
2683 0 : h2[0] = GetLayer(iL2)->GetHit(iH2)->GetPos()[0];
2684 0 : h2[1] = GetLayer(iL2)->GetHit(iH2)->GetPos()[1];
2685 0 : h2[2] = GetLayer(iL2)->GetHit(iH2)->GetPos()[2];
2686 :
2687 0 : TVector3 vec(h2[0]-h1[0],h2[1]-h1[1],h2[2]-h1[2]);
2688 :
2689 0 : 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 0 : if (!RuleSelectCell(h1,h2,iL1)) continue;
2694 0 : 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 0 : for (Int_t i = 0; i < 5; i++) hitSta[i] = kFALSE;
2699 : nHitSta = 0;
2700 :
2701 0 : hitSta[iL1/2] = kTRUE;
2702 0 : hitSta[iL2/2] = kTRUE;
2703 :
2704 : nptr = 0;
2705 :
2706 0 : trX[nptr] = h1[0];
2707 0 : trY[nptr] = h1[1];
2708 0 : trZ[nptr] = h1[2];
2709 0 : lay[nptr] = iL1;
2710 0 : hit[nptr] = iH1;
2711 0 : trkid[nptr] = GetLayer(iL1)->GetHit(iH1)->GetTrackGID();
2712 0 : detid[nptr] = GetLayer(iL1)->GetHit(iH1)->GetDetElemID();
2713 : nptr++;
2714 :
2715 0 : for (Int_t iL = (iL1+1); iL <= (iL2-1); iL++) {
2716 :
2717 0 : nH = GetLayer(iL)->GetNhits();
2718 :
2719 0 : if (prn) printf("L %d nH %d \n",iL,nH);
2720 :
2721 : dRmin = dRcut;
2722 : addHit = kFALSE;
2723 :
2724 0 : for (Int_t iH = 0; iH < nH; iH++) {
2725 :
2726 : //if (prn) printf("H %d \n",iH);
2727 :
2728 0 : if (GetLayer(iL)->GetHit(iH)->IsUsed()) continue;
2729 :
2730 0 : h[0] = GetLayer(iL)->GetHit(iH)->GetPos()[0];
2731 0 : h[1] = GetLayer(iL)->GetHit(iH)->GetPos()[1];
2732 0 : h[2] = GetLayer(iL)->GetHit(iH)->GetPos()[2];
2733 :
2734 0 : hx = h1[0] + (h2[0]-h1[0])*(h[2]-h1[2])/(h2[2]-h1[2]);
2735 0 : hy = h1[1] + (h2[1]-h1[1])*(h[2]-h1[2])/(h2[2]-h1[2]);
2736 :
2737 0 : dR = TMath::Sqrt((hx-h[0])*(hx-h[0])+(hy-h[1])*(hy-h[1]));
2738 0 : if (dR >= dRmin) continue;
2739 0 : AliDebug(1,Form("Hit %d added",iH));
2740 0 : hitSta[iL/2] = kTRUE;
2741 :
2742 : dRmin = dR;
2743 :
2744 0 : trX[nptr] = h[0];
2745 0 : trY[nptr] = h[1];
2746 0 : trZ[nptr] = h[2];
2747 0 : lay[nptr] = iL;
2748 0 : hit[nptr] = iH;
2749 0 : trkid[nptr] = GetLayer(iL)->GetHit(iH)->GetTrackGID();
2750 0 : detid[nptr] = GetLayer(iL)->GetHit(iH)->GetDetElemID();
2751 :
2752 : addHit = kTRUE;
2753 :
2754 0 : } // loop hits intermediate layer
2755 :
2756 0 : if (addHit) nptr++;
2757 :
2758 : } // loop intermediate layers
2759 :
2760 0 : trX[nptr] = h2[0];
2761 0 : trY[nptr] = h2[1];
2762 0 : trZ[nptr] = h2[2];
2763 0 : lay[nptr] = iL2;
2764 0 : hit[nptr] = iH2;
2765 0 : trkid[nptr] = GetLayer(iL2)->GetHit(iH2)->GetTrackGID();
2766 0 : detid[nptr] = GetLayer(iL2)->GetHit(iH2)->GetDetElemID();
2767 0 : nptr++;
2768 :
2769 : // create cells and tracks
2770 : //
2771 : // calculate how many stations have hit(s)
2772 0 : AliDebug(2,Form("nptr %d fMinTrackLength %d",nptr,fMinTrackLength));
2773 0 : if (nptr < fMinTrackLength) continue;
2774 : nHitSta = 0;
2775 0 : for (Int_t i = 0; i < 5; i++) {
2776 0 : if (hitSta[i]) nHitSta++;
2777 : }
2778 0 : AliDebug(2,Form("nHitSta %d fMinTrackLength %d",nHitSta,fMinTrackLength));
2779 0 : 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 0 : AliDebug(1,Form("Adding Track %d ",trackGID));
2793 0 : track = AddTrack(trackGID++);
2794 0 : track->SetStartLayer(iL1);
2795 :
2796 : // loop over hits in reverse order, like in RunBackward
2797 0 : for (Int_t iptr = (nptr-1); iptr >= 1; iptr--) {
2798 :
2799 0 : trkid1 = trkid[iptr-1];
2800 0 : lay1 = lay[iptr-1];
2801 0 : hit1 = hit[iptr-1];
2802 0 : detid1 = detid[iptr-1];
2803 :
2804 0 : htr1[0] = GetLayer(lay1)->GetHit(hit1)->GetPos()[0];
2805 0 : htr1[1] = GetLayer(lay1)->GetHit(hit1)->GetPos()[1];
2806 0 : htr1[2] = GetLayer(lay1)->GetHit(hit1)->GetPos()[2];
2807 0 : mftClsId1 = GetLayer(lay1)->GetHit(hit1)->GetMFTClsId();
2808 :
2809 0 : GetLayer(lay1)->GetHit(hit1)->SetUsed();
2810 :
2811 0 : trkid2 = trkid[iptr];
2812 0 : lay2 = lay[iptr];
2813 0 : hit2 = hit[iptr];
2814 0 : detid2 = detid[iptr];
2815 :
2816 0 : htr2[0] = GetLayer(lay2)->GetHit(hit2)->GetPos()[0];
2817 0 : htr2[1] = GetLayer(lay2)->GetHit(hit2)->GetPos()[1];
2818 0 : htr2[2] = GetLayer(lay2)->GetHit(hit2)->GetPos()[2];
2819 0 : mftClsId2 = GetLayer(lay2)->GetHit(hit2)->GetMFTClsId();
2820 :
2821 0 : GetLayer(lay2)->GetHit(hit2)->SetUsed();
2822 :
2823 : // create cells
2824 0 : cell = GetLayer(lay[iptr-1])->AddCell();
2825 0 : cell->SetHits(htr1,htr2,fPlanesZ[lay1],fPlanesZ[lay2]);
2826 0 : cell->SetMFTClsId(mftClsId1,mftClsId2);
2827 0 : cell->SetLayers(lay1,lay2);
2828 0 : cell->SetStatus(1);
2829 0 : cell->SetGID(fCellGID++,trkid1,trkid2);
2830 0 : cell->SetDetElemID(detid1,detid2);
2831 :
2832 : // add cell to track
2833 0 : track->AddCell(cell);
2834 0 : cell->SetUsed(kTRUE);
2835 :
2836 : }
2837 :
2838 0 : } // loop hits last layer
2839 0 : } // loop hits first layer
2840 :
2841 : } // end seed
2842 :
2843 0 : }
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 0 : Double_t xCl[nMaxh], yCl[nMaxh], yErr[nMaxh];
2859 : Int_t idet = 0;
2860 0 : for (Int_t i = 0; i < nDet; i++) {
2861 0 : if (i == skip) continue;
2862 0 : xCl[idet] = xcl[i];
2863 0 : yCl[idet] = ycl[i];
2864 0 : yErr[idet] = yerr[i];
2865 0 : idet++;
2866 0 : }
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 0 : for (Int_t i = 0; i < idet; i++) {
2874 0 : S1 += 1.0/(yErr[i]*yErr[i]);
2875 0 : SXY += xCl[i]*yCl[i]/(yErr[i]*yErr[i]);
2876 0 : SX += xCl[i]/(yErr[i]*yErr[i]);
2877 0 : SY += yCl[i]/(yErr[i]*yErr[i]);
2878 0 : SXX += xCl[i]*xCl[i]/(yErr[i]*yErr[i]);
2879 0 : if (i > 0) difx += TMath::Abs(xCl[i]-xCl[i-1]);
2880 0 : Xm += xCl[i];
2881 0 : Ym += yCl[i];
2882 0 : SsXX += xCl[i]*xCl[i];
2883 0 : SsYY += yCl[i]*yCl[i];
2884 0 : SsXY += xCl[i]*yCl[i];
2885 : }
2886 0 : delta = SXX*S1 - SX*SX;
2887 0 : if (delta == 0.) {
2888 0 : return kFALSE;
2889 : }
2890 0 : a = (SXY*S1 - SX*SY)/delta;
2891 0 : b = (SY*SXX - SX*SXY)/delta;
2892 :
2893 0 : Ym /= (Double_t)idet;
2894 0 : Xm /= (Double_t)idet;
2895 0 : SsYY -= (Double_t)idet*(Ym*Ym);
2896 0 : SsXX -= (Double_t)idet*(Xm*Xm);
2897 0 : SsXY -= (Double_t)idet*(Ym*Xm);
2898 : Double_t eps = 1.E-24;
2899 0 : if ((idet > 2) && (TMath::Abs(difx) > eps) && ((SsYY-(SsXY*SsXY)/SsXX) > 0.)) {
2900 0 : s = TMath::Sqrt((SsYY-(SsXY*SsXY)/SsXX)/(idet-2));
2901 0 : be = s*TMath::Sqrt(1./(Double_t)idet+(Xm*Xm)/SsXX);
2902 0 : ae = s/TMath::Sqrt(SsXX);
2903 0 : } else {
2904 0 : be = 0.;
2905 0 : ae = 0.;
2906 : }
2907 :
2908 : return kTRUE;
2909 :
2910 0 : }
2911 :
|
