]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - ITS/AliITStrackerSA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
typo fix
[u/mrichter/AliRoot.git] / ITS / AliITStrackerSA.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-2003, 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/* $Id$ */
17
18///////////////////////////////////////////////////////////
19// Stand alone ITS tracker class //
20// Origin: Elisabetta Crescio - crescio@to.infn.it //
21// Updated: Francesco Prino - prino@to.infn.it //
22///////////////////////////////////////////////////////////
23
24#include <stdlib.h>
25
26#include <TArrayI.h>
27#include <TBranch.h>
28#include <TObjArray.h>
29#include <TTree.h>
30
31#include "AliESDEvent.h"
32#include "AliESDVertex.h"
33#include "AliESDtrack.h"
34#include "AliITSVertexer.h"
35#include "AliITSclusterTable.h"
36#include "AliITSRecPoint.h"
37#include "AliITSgeomTGeo.h"
38#include "AliITStrackSA.h"
39#include "AliITStrackerSA.h"
40#include "AliITSReconstructor.h"
41#include "AliLog.h"
42#include "AliRun.h"
43
44ClassImp(AliITStrackerSA)
45
46//____________________________________________________________________________
47AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
48fPhiEstimate(0),
49fITSStandAlone(0),
50fLambdac(0),
51fPhic(0),
52fCoef1(0),
53fCoef2(0),
54fCoef3(0),
55fNloop(0),
56fPhiWin(0),
57fLambdaWin(0),
58fListOfTracks(0),
59fListOfSATracks(0),
60fITSclusters(0),
61fInwardFlag(0),
62fOuterStartLayer(0),
63fInnerStartLayer(5),
64fMinNPoints(0),
65fMinQ(0.),
66fCluCoord(0){
67 // Default constructor
68 Init();
69
70}
71//____________________________________________________________________________
72AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
73fPhiEstimate(0),
74fITSStandAlone(0),
75fLambdac(0),
76fPhic(0),
77fCoef1(0),
78fCoef2(0),
79fCoef3(0),
80fNloop(0),
81fPhiWin(0),
82fLambdaWin(0),
83fListOfTracks(0),
84fListOfSATracks(0),
85fITSclusters(0),
86fInwardFlag(0),
87fOuterStartLayer(0),
88fInnerStartLayer(5),
89fMinNPoints(0),
90fMinQ(0.),
91fCluCoord(0)
92{
93 // Standard constructor (Vertex is known and passed to this obj.)
94 if (geom) {
95 AliWarning("\"geom\" is actually a dummy argument !");
96 }
97
98 Init();
99
100}
101
102//____________________________________________________________________________
103AliITStrackerSA::~AliITStrackerSA(){
104 // destructor
105
106 if(fPhiWin)delete []fPhiWin;
107 if(fLambdaWin)delete []fLambdaWin;
108 fListOfTracks->Delete();
109 delete fListOfTracks;
110 fListOfSATracks->Delete();
111 delete fListOfSATracks;
112 if(fCluCoord){
113 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
114 if(fCluCoord[i]){
115 fCluCoord[i]->Delete();
116 delete fCluCoord[i];
117 }
118 }
119 delete [] fCluCoord;
120 }
121}
122
123//____________________________________________________________________________
124Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
125// This method is used to find and fit the tracks. By default the corresponding
126// method in the parent class is invoked. In this way a combined tracking
127// TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
128// is done in the ITS only. In the standard reconstruction chain this option
129// can be set via AliReconstruction::SetOption("ITS","onlyITS")
130 Int_t rc=0;
131
132 if(!fITSStandAlone){
133 rc=AliITStrackerMI::Clusters2Tracks(event);
134 }
135 else {
136 AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
137 }
138 if(!rc){
139 rc=FindTracks(event,kFALSE);
140 Int_t nSPDcontr=0;
141 const AliESDVertex *spdv = event->GetPrimaryVertexSPD();
142 if(spdv) nSPDcontr = spdv->GetNContributors();
143 if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE &&
144 nSPDcontr<=AliITSReconstructor::GetRecoParam()->GetMaxSPDcontrForSAToUseAllClusters()) {
145 rc=FindTracks(event,kTRUE);
146 }
147 }
148 return rc;
149}
150
151//____________________________________________________________________________
152void AliITStrackerSA::Init(){
153 // Reset all data members
154 fPhiEstimate=0;
155 for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
156 fLambdac=0;
157 fPhic=0;
158 fCoef1=0;
159 fCoef2=0;
160 fCoef3=0;
161 fPointc[0]=0;
162 fPointc[1]=0;
163 Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
164 if(nLoops==32){
165 SetFixedWindowSizes();
166 }else{
167 Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
168 Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
169 Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
170 Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
171 SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
172 }
173 fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
174 fITSclusters = 0;
175 SetOuterStartLayer(1);
176 SetSAFlag(kFALSE);
177 fListOfTracks=new TClonesArray("AliITStrackMI",100);
178 fListOfSATracks=new TClonesArray("AliITStrackSA",100);
179 fCluCoord = 0;
180 fMinNPoints = 3;
181 }
182//_______________________________________________________________________
183void AliITStrackerSA::ResetForFinding(){
184 // Reset data members used in all loops during track finding
185 fPhiEstimate=0;
186 for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
187 fLambdac=0;
188 fPhic=0;
189 fCoef1=0;
190 fCoef2=0;
191 fCoef3=0;
192 fPointc[0]=0;
193 fPointc[1]=0;
194 fListOfTracks->Clear();
195 fListOfSATracks->Clear();
196}
197
198
199
200//______________________________________________________________________
201Int_t AliITStrackerSA::FindTracks(AliESDEvent* event, Bool_t useAllClusters){
202
203// Track finder using the ESD object
204
205 AliDebug(2,Form(" field is %f",event->GetMagneticField()));
206 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
207
208 if(!fITSclusters){
209 Fatal("FindTracks","ITS cluster tree is not accessed!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
210 return -1;
211 }
212 //Reads event and mark clusters of traks already found, with flag kITSin
213 Int_t nentr=event->GetNumberOfTracks();
214 if(!useAllClusters) {
215 while (nentr--) {
216 AliESDtrack *track=event->GetTrack(nentr);
217 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
218 Int_t idx[12];
219 Int_t ncl = track->GetITSclusters(idx);
220 for(Int_t k=0;k<ncl;k++){
221 AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
222 cll->SetBit(kSAflag);
223 }
224 }
225 }
226 }else{
227 while (nentr--) {
228 AliESDtrack *track=event->GetTrack(nentr);
229 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
230 Int_t idx[12];
231 Int_t ncl = track->GetITSclusters(idx);
232 for(Int_t k=0;k<ncl;k++){
233 AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
234 cll->ResetBit(kSAflag);
235 }
236 }
237 }
238 }
239 //Get primary vertex
240 Double_t primaryVertex[3];
241 event->GetVertex()->GetXYZ(primaryVertex);
242 //Creates TClonesArray with clusters for each layer. The clusters already used
243 //by AliITStrackerMI are not considered
244 Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
245 Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
246 if (fCluCoord == 0) {
247 fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
248 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
249 fCluCoord[i]=0;
250 }
251 }
252 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
253 AliITSlayer &layer=fgLayers[i];
254 if (!ForceSkippingOfLayer(i)) {
255 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
256 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
257 if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
258 if(cls->GetQ()==0) continue; //fake clusters dead zones
259 if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
260 nclusters[i]++;
261 }
262 }
263 dmar[i]=0;
264 if(!fCluCoord[i]){
265 fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
266 }else{
267 fCluCoord[i]->Delete();
268 fCluCoord[i]->Expand(nclusters[i]);
269 }
270 }
271
272 for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
273 TClonesArray &clucoo = *fCluCoord[ilay];
274 AliITSlayer &layer=fgLayers[ilay];
275 if (!ForceSkippingOfLayer(ilay)) {
276 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
277 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
278 if(cls->TestBit(kSAflag)==kTRUE) continue;
279 if(cls->GetQ()==0) continue;
280 if(ilay>1 && cls->GetQ()<=fMinQ) continue;
281 Double_t phi=0;Double_t lambda=0;
282 Double_t x=0;Double_t y=0;Double_t z=0;
283 Double_t sx=0;Double_t sy=0;Double_t sz=0;
284 GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
285 GetCoorErrors(cls,sx,sy,sz);
286 new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
287 dmar[ilay]++;
288 }
289 }
290 fCluCoord[ilay]->Sort();
291 }
292
293 // track counter
294 Int_t ntrack=0;
295
296 static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
297 Int_t startLayForSeed=0;
298 Int_t lastLayForSeed=fOuterStartLayer;
299 Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
300 Int_t seedStep=1;
301 if(fInwardFlag){
302 startLayForSeed=AliITSgeomTGeo::GetNLayers()-1;
303 lastLayForSeed=fInnerStartLayer;
304 nSeedSteps=startLayForSeed-lastLayForSeed;
305 seedStep=-1;
306 }
307
308 // loop on minimum number of points
309 for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) {
310
311 // loop on starting layer for track finding
312 for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
313 Int_t theLay=startLayForSeed+iSeedLay*seedStep;
314 if(ForceSkippingOfLayer(theLay)) continue;
315 Int_t minNPoints=iMinNPoints-theLay;
316 if(fInwardFlag) minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-theLay);
317 for(Int_t i=theLay+1;i<AliITSgeomTGeo::GetNLayers();i++)
318 if(ForceSkippingOfLayer(i))
319 minNPoints--;
320 if(minNPoints<fMinNPoints) continue;
321
322 // loop on phi and lambda window size
323 for(Int_t nloop=0;nloop<fNloop;nloop++){
324 Int_t nclTheLay=fCluCoord[theLay]->GetEntries();
325 while(nclTheLay--){
326 ResetForFinding();
327 Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
328 if(!useRP) continue;
329 AliITStrackSA trs;
330
331 Int_t pflag=0;
332 Int_t kk;
333 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
334
335 kk=0;
336 nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
337 &trs,primaryVertex[2],pflag);
338 Int_t nextLay=theLay+seedStep;
339 Bool_t goon=kTRUE;
340 if(nextLay<0 || nextLay == 6) goon = kFALSE;
341 while(goon){
342 kk++;
343 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
344 &trs,primaryVertex[2],pflag);
345 if(nClusLay[kk]!=0){
346 pflag=1;
347 if(kk==1) {
348 fPoint3[0]=fPointc[0];
349 fPoint3[1]=fPointc[1];
350 } else {
351 UpdatePoints();
352 }
353 }
354 nextLay+=seedStep;
355 if(nextLay<0 || nextLay==6) goon=kFALSE;
356 }
357
358
359 Int_t layOK=0;
360 if(!fInwardFlag){
361 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-theLay;nnp++){
362 if(nClusLay[nnp]!=0) layOK+=1;
363 }
364 }else{
365 for(Int_t nnp=theLay; nnp>=0; nnp--){
366 if(nClusLay[nnp]!=0) layOK+=1;
367 }
368 }
369 if(layOK>=minNPoints){
370 AliDebug(2,Form("---NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]));
371 AliITStrackV2* tr2 = 0;
372 tr2 = FitTrack(&trs,primaryVertex);
373 if(!tr2){
374 continue;
375 }
376 AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
377
378 StoreTrack(tr2,event,useAllClusters);
379 ntrack++;
380
381 }
382
383 }//end loop on clusters of theLay
384 } //end loop on window sizes
385 } //end loop on theLay
386 }//end loop on min points
387
388 // search for 1-point tracks in SPD, only for cosmics
389 // (A.Dainese 21.03.08)
390 if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
391 TMath::Abs(event->GetMagneticField())<0.01) {
392 Int_t outerLayer=1; // only SPD
393 for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
394 // counter for clusters on each layer
395
396 for(Int_t nloop=0;nloop<fNloop;nloop++){
397 Int_t nclInnLay=fCluCoord[innLay]->GetEntries();
398 while(nclInnLay--){ //loop starting from layer innLay
399 ResetForFinding();
400 Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
401 if(!useRP) continue;
402 AliITStrackSA trs;
403
404 Int_t pflag=0;
405 Int_t kk;
406 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
407
408 kk=0;
409 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
410 &trs,primaryVertex[2],pflag);
411 for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
412 kk++;
413 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
414 &trs,primaryVertex[2],pflag);
415 if(nClusLay[kk]!=0){
416 pflag=1;
417 if(kk==1) {
418 fPoint3[0]=fPointc[0];
419 fPoint3[1]=fPointc[1];
420 } else {
421 UpdatePoints();
422 }
423 }
424 }
425
426 Int_t layOK=0;
427 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
428 if(nClusLay[nnp]!=0) layOK+=1;
429 }
430 if(layOK==1) {
431 AliDebug(2,Form("----NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]));
432 AliITStrackV2* tr2 = 0;
433 Bool_t onePoint = kTRUE;
434 tr2 = FitTrack(&trs,primaryVertex,onePoint);
435 if(!tr2){
436 continue;
437 }
438 AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
439
440 StoreTrack(tr2,event,useAllClusters);
441 ntrack++;
442
443 }
444
445 }//end loop on clusters of innLay
446 } //end loop on window sizes
447
448 } //end loop on innLay
449 } // end search 1-point tracks
450
451 if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
452 ResetForFinding();
453 return 0;
454
455}
456
457//________________________________________________________________________
458
459AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
460 //fit of the found track (most general case, also <6 points, layers missing)
461 // A.Dainese 16.11.07
462
463
464 const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
465
466 static Int_t firstmod[AliITSgeomTGeo::kNLayers];
467 static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
468 static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
469 static Int_t end[AliITSgeomTGeo::kNLayers];
470 static Int_t indices[AliITSgeomTGeo::kNLayers];
471
472 static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
473
474 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
475 firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
476 end[i]=0;
477 for(Int_t j=0;j<kMaxClu; j++){
478 clind[i][j]=0;
479 clmark[i][j]=0;
480 listlayer[i][j]=0;
481 }
482 }
483
484
485 Int_t nclusters = tr->GetNumberOfClustersSA();
486 for(Int_t ncl=0;ncl<nclusters;ncl++){
487 Int_t index = tr->GetClusterIndexSA(ncl);
488 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
489 Int_t lay = (index & 0xf0000000) >> 28;
490 Int_t nInLay=end[lay];
491 listlayer[lay][nInLay]=cl;
492 clind[lay][nInLay]=index;
493 end[lay]++;
494 }
495
496 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
497 for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
498 Int_t mark = tr->GetClusterMark(nlay,ncl);
499 clmark[nlay][ncl]=mark;
500 }
501 }
502
503
504 Int_t firstLay=-1,secondLay=-1;
505 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
506 if(end[i]==0) {
507 end[i]=1;
508 }else{
509 if(firstLay==-1) {
510 firstLay=i;
511 } else if(secondLay==-1) {
512 secondLay=i;
513 }
514 }
515 }
516
517 if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
518 TClonesArray &arrMI= *fListOfTracks;
519 TClonesArray &arrSA= *fListOfSATracks;
520 Int_t nFoundTracks=0;
521
522
523 for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
524 indices[0]=l0;
525 for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
526 indices[1]=l1;
527 for(Int_t l2=0;l2<end[2];l2++){ //loop on layer 3
528 indices[2]=l2;
529 for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
530 indices[3]=l3;
531 for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
532 indices[4]=l4;
533 for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
534 indices[5]=l5;
535
536 // estimate curvature from 2 innermost points (or innermost point + vertex)
537
538 Int_t iFirstLay=indices[firstLay];
539 Int_t mrk1=clmark[firstLay][iFirstLay];
540
541 AliITSRecPoint* p1=(AliITSRecPoint*)listlayer[firstLay][iFirstLay];
542 Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
543 Int_t layer,ladder,detector;
544 AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
545 Double_t yclu1 = p1->GetY();
546 Double_t zclu1 = p1->GetZ();
547
548 Double_t x1,y1,z1;
549 Double_t x2,y2,z2;
550 Double_t cv=0,tgl2=0,phi2=0;
551 AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
552 x1 = arr1->GetX();
553 y1 = arr1->GetY();
554 z1 = arr1->GetZ();
555
556 if(secondLay>0) {
557 Int_t iSecondLay=indices[secondLay];
558 Int_t mrk2=clmark[secondLay][iSecondLay];
559 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
560 x2 = arr2->GetX();
561 y2 = arr2->GetY();
562 z2 = arr2->GetZ();
563 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
564 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
565 phi2 = TMath::ATan2((y2-y1),(x2-x1));
566 } else { // special case of 1-point tracks, only for cosmics (B=0)
567 x2 = primaryVertex[0];
568 y2 = primaryVertex[1];
569 z2 = primaryVertex[2];
570 cv = 0;
571 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
572 phi2 = TMath::ATan2((y1-y2),(x1-x2));
573 }
574
575 // create track and attach it the RecPoints
576 AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
577 for(Int_t iLay=5; iLay>=0; iLay--){
578 Int_t iInLay=indices[iLay];
579 AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay];
580 if(cl!=0){
581 trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
582 trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
583 }
584 }
585
586 //fit with Kalman filter using AliITStrackerMI::RefitAt()
587 AliITStrackSA ot(trac);
588
589 ot.ResetCovariance(10.);
590 ot.ResetClusters();
591
592 // Propagate inside the innermost layer with a cluster
593 if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
594
595 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6
596 AliITStrackMI otrack2(ot);
597 otrack2.ResetCovariance(10.);
598 otrack2.ResetClusters();
599 //fit from layer 6 to layer 1
600 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) {
601 new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
602 new(arrSA[nFoundTracks]) AliITStrackSA(trac);
603 ++nFoundTracks;
604 }
605
606 }
607 }
608 }//end loop layer 6
609 }//end loop layer 5
610 }//end loop layer 4
611 }//end loop layer 3
612 }//end loop layer 2
613 }//end loop layer 1
614
615
616
617
618 if(fListOfTracks->GetEntries()==0) return 0;
619
620 Int_t lowchi2 = FindTrackLowChiSquare();
621 AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
622 AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
623
624 if(otrack==0) return 0;
625
626 CookLabel(otrack,0.); //MI change - to see fake ratio
627 Int_t label=FindLabel(otrack);
628 otrack->SetLabel(label);
629 Double_t low=0.;
630 Double_t up=0.51;
631 otrack->CookdEdx(low,up);
632
633 //remove clusters of found track
634 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
635 for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
636 Int_t index = trsa->GetClusterMark(nlay,cln);
637 RemoveClusterCoord(nlay,index);
638 }
639 }
640
641 return otrack;
642
643}
644
645//_______________________________________________________
646void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const
647{
648 //
649 // Add new track to the ESD
650 //
651 AliESDtrack outtrack;
652 outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
653 if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
654 for(Int_t i=0;i<12;i++) {
655 outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
656 }
657 Double_t sdedx[4]={0.,0.,0.,0.};
658 for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
659 outtrack.SetITSdEdxSamples(sdedx);
660
661
662 if(AliITSReconstructor::GetRecoParam()->GetSAUsedEdxInfo()){
663 Double_t mom=t->P();
664 Double_t ppid[AliPID::kSPECIES];
665 for(Int_t isp=0;isp<AliPID::kSPECIES;isp++) ppid[isp]=0.;
666 ppid[AliPID::kPion]=1.;
667 if(mom<0.7){
668 Double_t truncmean=t->GetdEdx();
669 Int_t ide=fITSPid->GetParticleIdFromdEdxVsP(mom,truncmean,kTRUE);
670 if(ide==AliPID::kProton){
671 ppid[AliPID::kProton]=1.;
672 ppid[AliPID::kPion]=0.;
673 }
674 else if(ide==AliPID::kKaon){
675 ppid[AliPID::kKaon]=1.;
676 ppid[AliPID::kPion]=0.;
677 }
678 }
679 outtrack.SetITSpid(ppid);
680 outtrack.SetESDpid(ppid);
681 }
682 event->AddTrack(&outtrack);
683
684 return;
685}
686
687
688//_______________________________________________________
689Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
690 //function used to to find the clusters associated to the track
691
692 if(ForceSkippingOfLayer(layer)) return 0;
693
694
695 Int_t nc=0;
696 AliITSlayer &lay = fgLayers[layer];
697 Double_t r=lay.GetR();
698 if(pflag==1){
699 Double_t cx1,cx2,cy1,cy2;
700 FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
701 if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
702 return 0;
703 Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
704 Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
705 fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
706 }
707
708
709 Double_t phiExpect=fPhiEstimate;
710 Double_t lamExpect=fLambdac;
711
712 Int_t ncl = fCluCoord[layer]->GetEntriesFast();
713 Int_t startcl=FindIndex(layer,lamExpect-lambdawindow*1.02);
714 Int_t endcl=FindIndex(layer,lamExpect+lambdawindow*1.02)+1;
715 if(endcl>=ncl) endcl=ncl-1;
716
717 for (Int_t index=startcl; index<=endcl; index++) {
718 //for (Int_t index=0; index<ncl; index++) {
719 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
720
721
722 Double_t phi = arr->GetPhi();
723 Double_t deltaPhi = phi-phiExpect;
724 if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
725 else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
726 if (TMath::Abs(deltaPhi)>phiwindow) continue;
727
728 Double_t lambda = arr->GetLambda();
729 if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
730
731 if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
732 if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
733 Int_t orind = arr->GetOrInd();
734 trs->AddClusterSA(layer,orind);
735 trs->AddClusterMark(layer,index);
736 nc++;
737 fLambdac=lambda;
738 fPhiEstimate=phi;
739
740 fPointc[0]=arr->GetX();
741 fPointc[1]=arr->GetY();
742
743 }
744 return nc;
745}
746
747//________________________________________________________________
748Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
749 // Sets the first point (seed) for tracking
750
751 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
752 fPhic = arr->GetPhi();
753 fLambdac = arr->GetLambda();
754 fPhiEstimate = fPhic;
755 fPoint1[0]=primaryVertex[0];
756 fPoint1[1]=primaryVertex[1];
757 fPoint2[0]=arr->GetX();
758 fPoint2[1]=arr->GetY();
759 return kTRUE;
760}
761
762//________________________________________________________________
763void AliITStrackerSA::UpdatePoints(){
764 //update of points for the estimation of the curvature
765
766 fPoint2[0]=fPoint3[0];
767 fPoint2[1]=fPoint3[1];
768 fPoint3[0]=fPointc[0];
769 fPoint3[1]=fPointc[1];
770
771
772}
773
774//___________________________________________________________________
775Int_t AliITStrackerSA::FindEquation(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3,Double_t& a, Double_t& b, Double_t& c){
776
777 //given (x,y) of three recpoints (in global coordinates)
778 //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
779 double dx31=x3-x1,dy31=y3-y1,dx21=x2-x1,dy21=y2-y1;
780 Double_t den = dx31*dy21-dx21*dy31;
781 if(den==0) return 0;
782 else den = 1./den;
783 //
784 double r31 = -dx31*(x1+x3) - dy31*(y1+y3);
785 double r21 = -dx21*(x1+x2) - dy21*(y1+y2);
786 //
787 double da = r31*dy21 - r21*dy31;
788 double db = r21*dx31 - r31*dx21;
789 a = da*den;
790 b = db*den;
791 c = -x1*x1-y1*y1-a*x1-b*y1;
792 return 1;
793 }
794//__________________________________________________________________________
795 Int_t AliITStrackerSA::FindIntersection(Double_t a1, Double_t b1, Double_t c1, Double_t c2,Double_t& x1,Double_t& y1, Double_t& x2, Double_t& y2){
796
797 //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
798 //c2 is -rlayer*rlayer
799
800 if(a1==0) return 0;
801 Double_t m = c2-c1;
802 Double_t aA = (b1*b1)/(a1*a1)+1;
803 Double_t bB = (-2*m*b1/(a1*a1));
804 Double_t cC = c2+(m*m)/(a1*a1);
805 Double_t dD = bB*bB-4*aA*cC;
806 if(dD<0) return 0;
807
808 y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
809 y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
810 x1 = (c2-c1-b1*y1)/a1;
811 x2 = (c2-c1-b1*y2)/a1;
812
813 return 1;
814}
815//____________________________________________________________________
816Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
817x2,Double_t y2,Double_t x3,Double_t y3){
818
819 //calculates the curvature of track
820 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
821 if(den==0) return 0;
822 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
823 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
824 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
825 Double_t xc=-a/2.;
826
827 if((a*a+b*b-4*c)<0) return 0;
828 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
829 if(rad==0) return 0;
830
831 if((x1>0 && y1>0 && x1<xc)) rad*=-1;
832 if((x1<0 && y1>0 && x1<xc)) rad*=-1;
833 // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
834 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
835
836 return 1/rad;
837
838}
839
840
841//____________________________________________________________________
842Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
843
844 //Returns the point closest to pp
845
846 Double_t diff1 = p1-pp;
847 Double_t diff2 = p2-pp;
848
849 if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
850 else fPhiEstimate=p2;
851 return fPhiEstimate;
852
853}
854
855
856//_________________________________________________________________
857Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
858 // returns track with lowest chi square
859 Int_t dim=fListOfTracks->GetEntries();
860 if(dim<=1) return 0;
861 AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
862 Double_t minChi2=trk->GetChi2();
863 Int_t index=0;
864 for(Int_t i=1;i<dim;i++){
865 trk = (AliITStrackV2*)fListOfTracks->At(i);
866 Double_t chi2=trk->GetChi2();
867 if(chi2<minChi2){
868 minChi2=chi2;
869 index=i;
870 }
871 }
872 return index;
873}
874
875//__________________________________________________________
876Int_t AliITStrackerSA::FindLabel(AliITStrackV2* track){
877 // compute the track label starting from cluster labels
878
879 Int_t labl[AliITSgeomTGeo::kNLayers][3];
880 Int_t cnts[AliITSgeomTGeo::kNLayers][3];
881 for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
882 for(Int_t k=0;k<3;k++){
883 labl[j][k]=-2;
884 cnts[j][k]=1;
885 }
886 }
887 Int_t iNotLabel=0;
888 for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
889 Int_t indexc = track->GetClusterIndex(i);
890 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
891 Int_t iLayer=cl->GetLayer();
892 for(Int_t k=0;k<3;k++){
893 labl[iLayer][k]=cl->GetLabel(k);
894 if(labl[iLayer][k]<0) iNotLabel++;
895 }
896 }
897 if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
898
899 for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
900 for(Int_t j2=0; j2<j1; j2++){
901 for(Int_t k1=0; k1<3; k1++){
902 for(Int_t k2=0; k2<3; k2++){
903 if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
904 cnts[j2][k2]++;
905 cnts[j1][k1]=0;
906 }
907 }
908 }
909 }
910 }
911
912
913 Int_t cntMax=0;
914 Int_t label=-1;
915 for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
916 for(Int_t k=0;k<3;k++){
917 if(cnts[j][k]>cntMax && labl[j][k]>=0){
918 cntMax=cnts[j][k];
919 label=labl[j][k];
920 }
921 }
922 }
923
924 Int_t lflag=0;
925 for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
926 if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
927
928 if(lflag<track->GetNumberOfClusters()) label = -label;
929 return label;
930}
931//_____________________________________________________________________________
932void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Double_t phimin, Double_t phimax, Double_t lambdamin, Double_t lambdamax){
933 // Set sizes of the phi and lambda windows used for track finding
934 fNloop = n;
935 if(fPhiWin) delete [] fPhiWin;
936 if(fLambdaWin) delete [] fLambdaWin;
937 fPhiWin = new Double_t[fNloop];
938 fLambdaWin = new Double_t[fNloop];
939 Double_t stepPhi=(phimax-phimin)/(Double_t)(fNloop-1);
940 Double_t stepLambda=(lambdamax-lambdamin)/(Double_t)(fNloop-1);
941 for(Int_t k=0;k<fNloop;k++){
942 Double_t phi=phimin+k*stepPhi;
943 Double_t lam=lambdamin+k*stepLambda;
944 fPhiWin[k]=phi;
945 fLambdaWin[k]=lam;
946 }
947}
948//_____________________________________________________________________________
949void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
950 // Set sizes of the phi and lambda windows used for track finding
951 fNloop = n;
952 if(phi){ // user defined values
953 fPhiWin = new Double_t[fNloop];
954 fLambdaWin = new Double_t[fNloop];
955 for(Int_t k=0;k<fNloop;k++){
956 fPhiWin[k]=phi[k];
957 fLambdaWin[k]=lam[k];
958 }
959 }
960 else { // default values
961
962 Double_t phid[32] = {0.002,0.003,0.004,0.0045,0.0047,
963 0.005,0.0053,0.0055,0.006,0.0063,
964 0.0065,0.007,0.0073,0.0075,0.0077,
965 0.008,0.0083,0.0085,0.0087,0.009,
966 0.0095,0.0097,0.01,0.0105,0.011,
967 0.0115,0.012,0.0125,0.013,0.0135,
968 0.0140,0.0145};
969 Double_t lambdad[32] = {0.003,0.004,0.005,0.005,0.005,
970 0.005,0.005,0.006,0.006,0.006,
971 0.006,0.007,0.007,0.007,0.007,
972 0.007,0.007,0.007,0.007,0.007,
973 0.007,0.007,0.008,0.008,0.008,
974 0.008,0.008,0.008,0.008,0.008,
975 0.008,0.008};
976
977 if(fNloop!=32){
978 fNloop = 32;
979 }
980
981
982 fPhiWin = new Double_t[fNloop];
983 fLambdaWin = new Double_t[fNloop];
984
985 Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
986
987 for(Int_t k=0;k<fNloop;k++){
988 fPhiWin[k]=phid[k]*factor;
989 fLambdaWin[k]=lambdad[k]*factor;
990 }
991
992 }
993
994}
995//_______________________________________________________________________
996void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z, const Double_t* vertex){
997 //Returns values of phi (azimuthal) and lambda angles for a given cluster
998/*
999 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1000 Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1001 Double_t tx,ty,tz; fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1002
1003 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1004 Double_t phi1=TMath::Pi()/2+alpha;
1005 if (lay==1) phi1+=TMath::Pi();
1006
1007 Double_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1008 Double_t r=tx*cp+ty*sp;
1009
1010 xyz= r*cp - cl->GetY()*sp;
1011 y= r*sp + cl->GetY()*cp;
1012 z=cl->GetZ();
1013*/
1014 Float_t xyz[3];
1015 cl->GetGlobalXYZ(xyz);
1016 x=xyz[0];
1017 y=xyz[1];
1018 z=xyz[2];
1019
1020 phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1021 lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1022}
1023
1024//________________________________________________________________________
1025void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Double_t &sx,Double_t &sy, Double_t &sz){
1026
1027 //returns sigmax, y, z of cluster in global coordinates
1028/*
1029 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1030 Int_t lay,lad,det;
1031 AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1032
1033 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1034 Double_t phi=TMath::Pi()/2+alpha;
1035 if (lay==1) phi+=TMath::Pi();
1036
1037 Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1038*/
1039 Float_t covm[6];
1040 cl->GetGlobalCov(covm);
1041 sx=TMath::Sqrt(covm[0]);
1042 sy=TMath::Sqrt(covm[3]);
1043 sz=TMath::Sqrt(covm[5]);
1044/*
1045 sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1046 sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1047 sz = TMath::Sqrt(cl->GetSigmaZ2());
1048*/
1049}
1050
1051//________________________________________________________________________
1052Int_t AliITStrackerSA::FindIndex(Int_t lay, Double_t lamVal) const {
1053 // Find the cluster at limit of lambda window
1054
1055 Int_t base = 0;
1056 Int_t last = fCluCoord[lay]->GetEntriesFast()-1;
1057 if(last<0) return 0;
1058 Int_t position;
1059 Double_t lamfirst=((AliITSclusterTable*)fCluCoord[lay]->At(base))->GetLambda();
1060 if(lamfirst>lamVal) return base;
1061 Double_t lamlast=((AliITSclusterTable*)fCluCoord[lay]->At(last))->GetLambda();
1062 if(lamlast<=lamVal) return last;
1063 while (last >= base) {
1064 position = (base+last) / 2;
1065 Double_t a=((AliITSclusterTable*)fCluCoord[lay]->At(position))->GetLambda()-lamVal;
1066 Double_t b=((AliITSclusterTable*)fCluCoord[lay]->At(position+1))->GetLambda()-lamVal;
1067 if(a*b<=0) return position;
1068 if(a>0) last = position;
1069 else base = position;
1070 }
1071 return 0;
1072}
1073