1 /**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ////////////////////////////////////////////////////
19 // Stand alone tracker class //
20 // Origin: Elisabetta Crescio //
21 // e-mail: crescio@to.infn.it //
22 // tracks are saved as AliITStrackV2 objects //
23 ////////////////////////////////////////////////////
29 #include <TObjArray.h>
32 #include "AliESDEvent.h"
33 #include "AliESDVertex.h"
34 #include "AliESDtrack.h"
35 #include "AliITSVertexer.h"
36 #include "AliITSclusterTable.h"
37 #include "AliITSRecPoint.h"
38 #include "AliITSgeomTGeo.h"
39 #include "AliITStrackSA.h"
40 #include "AliITStrackerSA.h"
41 #include "AliITSReconstructor.h"
44 ClassImp(AliITStrackerSA)
46 //____________________________________________________________________________
47 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
68 // Default constructor
72 //____________________________________________________________________________
73 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
95 // Standard constructor (Vertex is known and passed to this obj.)
97 AliWarning("\"geom\" is actually a dummy argument !");
105 //____________________________________________________________________________
106 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliESDVertex *vert):AliITStrackerMI(0),
128 // Standard constructor (Vertex is known and passed to this obj.)
130 AliWarning("\"geom\" is actually a dummy argument !");
136 //____________________________________________________________________________
137 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliITSVertexer *vertexer):AliITStrackerMI(0),
159 // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
161 AliWarning("\"geom\" is actually a dummy argument !");
164 fVertexer = vertexer;
168 //____________________________________________________________________________
169 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(),
170 fPhiEstimate(tracker.fPhiEstimate),
171 fITSStandAlone(tracker.fITSStandAlone),
172 fLambdac(tracker.fLambdac),
173 fPhic(tracker.fPhic),
174 fCoef1(tracker.fCoef1),
175 fCoef2(tracker.fCoef2),
176 fCoef3(tracker.fCoef3),
177 fNloop(tracker.fNloop),
178 fPhiWin(tracker.fPhiWin),
179 fLambdaWin(tracker.fLambdaWin),
180 fVert(tracker.fVert),
181 fVertexer(tracker.fVertexer),
182 fListOfTracks(tracker.fListOfTracks),
183 fListOfSATracks(tracker.fListOfSATracks),
184 fITSclusters(tracker.fITSclusters),
185 fSixPoints(tracker.fSixPoints),
186 fOuterStartLayer(tracker.fOuterStartLayer),
187 fMinQ(tracker.fMinQ),
188 fCluLayer(tracker.fCluLayer),
189 fCluCoord(tracker.fCluCoord) {
191 for(Int_t i=0;i<2;i++){
192 fPoint1[i]=tracker.fPoint1[i];
193 fPoint2[i]=tracker.fPoint2[i];
194 fPoint3[i]=tracker.fPoint3[i];
195 fPointc[i]=tracker.fPointc[i];
197 if(tracker.fVertexer && tracker.fVert){
198 fVert = new AliESDVertex(*tracker.fVert);
201 fVert = tracker.fVert;
203 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
204 fCluLayer[i] = tracker.fCluLayer[i];
205 fCluCoord[i] = tracker.fCluCoord[i];
208 //______________________________________________________________________
209 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& source){
210 // Assignment operator.
211 this->~AliITStrackerSA();
212 new(this) AliITStrackerSA(source);
217 //____________________________________________________________________________
218 AliITStrackerSA::~AliITStrackerSA(){
220 // if fVertexer is not null, the AliESDVertex obj. is owned by this class
221 // and is deleted here
223 if(fVert)delete fVert;
228 if(fPhiWin)delete []fPhiWin;
229 if(fLambdaWin)delete []fLambdaWin;
230 fListOfTracks->Delete();
231 delete fListOfTracks;
232 fListOfSATracks->Delete();
233 delete fListOfSATracks;
235 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
237 fCluLayer[i]->Delete();
244 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
246 fCluCoord[i]->Delete();
255 //____________________________________________________________________________
256 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
257 // This method is used to find and fit the tracks. By default the corresponding
258 // method in the parent class is invoked. In this way a combined tracking
259 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
260 // is done in the ITS only. In the standard reconstruction chain this option
261 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
264 rc=AliITStrackerMI::Clusters2Tracks(event);
267 AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
269 if(!rc) rc=FindTracks(event);
273 //____________________________________________________________________________
274 void AliITStrackerSA::Init(){
275 // Reset all data members
277 for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
287 Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
289 SetFixedWindowSizes();
291 Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
292 Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
293 Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
294 Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
295 SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
297 fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
300 SetOuterStartLayer(0);
302 fListOfTracks=new TObjArray(0,0);
303 fListOfSATracks=new TObjArray(0,0);
307 //_______________________________________________________________________
308 void AliITStrackerSA::ResetForFinding(){
309 // Reset data members used in all loops during track finding
311 for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
319 fListOfTracks->Delete();
320 fListOfSATracks->Delete();
325 //______________________________________________________________________
326 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event){
328 // Track finder using the ESD object
333 Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
336 //Reads event and mark clusters of traks already found, with flag kITSin
337 Int_t nentr=event->GetNumberOfTracks();
338 if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kFALSE) {
340 AliESDtrack *track=event->GetTrack(nentr);
341 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
343 Int_t ncl = track->GetITSclusters(idx);
344 for(Int_t k=0;k<ncl;k++){
345 AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
346 cll->SetBit(kSAflag);
352 Double_t primaryVertex[3];
353 event->GetVertex()->GetXYZ(primaryVertex);
354 //Creates TClonesArray with clusters for each layer. The clusters already used
355 //by AliITStrackerMI are not considered
356 Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
357 Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
358 if (fCluLayer == 0) {
359 fCluLayer = new TClonesArray*[AliITSgeomTGeo::kNLayers];
360 fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
361 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
366 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
367 AliITSlayer &layer=fgLayers[i];
368 if (!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
369 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
370 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
371 if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
372 if(cls->GetQ()==0) continue; //fake clusters dead zones
373 if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
379 fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
381 fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
384 for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
385 TClonesArray &clulay = *fCluLayer[ilay];
386 TClonesArray &clucoo = *fCluCoord[ilay];
387 AliITSlayer &layer=fgLayers[ilay];
388 if (!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilay)) {
389 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
390 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
391 if(cls->TestBit(kSAflag)==kTRUE) continue;
392 if(cls->GetQ()==0) continue;
393 if(ilay>1 && cls->GetQ()<=fMinQ) continue;
394 Double_t phi=0;Double_t lambda=0;
395 Float_t x=0;Float_t y=0;Float_t z=0;
396 Float_t sx=0;Float_t sy=0;Float_t sz=0;
397 GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
398 GetCoorErrors(cls,sx,sy,sz);
399 new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
400 new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
409 Int_t minNPoints = (fSixPoints ? AliITSgeomTGeo::GetNLayers() : AliITSgeomTGeo::GetNLayers()-1);
410 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
411 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
416 static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
418 //loop on different minNPoints
419 Int_t minMinNPoints=minNPoints;
420 if(AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad()) minMinNPoints=2;
421 for(Int_t iMinNPoints=minNPoints; iMinNPoints>=minMinNPoints; iMinNPoints--) {
422 //loop on the different windows
423 for(Int_t nloop=0;nloop<fNloop;nloop++){
424 for(Int_t ncl=0;ncl<fCluLayer[0]->GetEntries();ncl++){ //loop starting from layer 0
429 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[0]->At(ncl);
432 if (cl->GetQ()<=0) continue;
434 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(0,ncl);
435 fPhic = arr->GetPhi();
436 fLambdac = arr->GetLambda();
437 fPhiEstimate = fPhic;
438 AliITStrackSA* trs = new AliITStrackSA();
439 fPoint1[0]=primaryVertex[0];
440 fPoint1[1]=primaryVertex[1];
443 fPoint2[0]=arr->GetX();
444 fPoint2[1]=arr->GetY();
445 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) nClusLay[i]=0;
446 nClusLay[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
447 nClusLay[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
450 fPoint3[0] = fPointc[0];
451 fPoint3[1] = fPointc[1];
453 nClusLay[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
454 if(nClusLay[1]==0 && nClusLay[2]==0) pflag=0;
455 if(nClusLay[2]!=0 && nClusLay[1]!=0){ pflag=1; UpdatePoints();}
456 if(nClusLay[2]!=0 && nClusLay[1]==0){
458 fPoint3[0]=fPointc[0];
459 fPoint3[1]=fPointc[1];
462 nClusLay[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
464 if(nClusLay[3]!=0) UpdatePoints();
465 nClusLay[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
467 if(nClusLay[4]!=0) UpdatePoints();
468 nClusLay[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
472 //check of the candidate track
473 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers();nnp++) {
474 if(nClusLay[nnp]!=0) layOK+=1;
477 if(layOK>=iMinNPoints) {
478 //printf("-NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
479 AliITStrackV2* tr2 = 0;
480 tr2 = FitTrack(trs,primaryVertex);
482 //printf("-NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
484 StoreTrack(tr2,event);
488 }//end loop on clusters of layer1
489 }//end loop on windows
490 }//end loop on min points
495 //if 5/6 points are required, second loop starting
496 //from second layer (SPD2), to find tracks with point of
497 //layer 1 missing. Not done for cosmics.
498 if(!fSixPoints && fOuterStartLayer==0) {
499 //printf("looking from SPD2\n");
500 // counter for clusters on each layer
501 for(Int_t nloop=0;nloop<fNloop;nloop++){
502 Int_t ncl2=fCluLayer[1]->GetEntries();
503 while(ncl2--){ //loop starting from layer 2
506 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[1]->At(ncl2);
509 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(1,ncl2);
510 fPhic = arr->GetPhi();
511 fLambdac = arr->GetLambda();
512 fPhiEstimate = fPhic;
514 AliITStrackSA* trs = new AliITStrackSA();
515 fPoint1[0]=primaryVertex[0];
516 fPoint1[1]=primaryVertex[1];
518 fPoint2[0]=arr->GetX();
519 fPoint2[1]=arr->GetY();
520 for(Int_t kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
521 nClusLay[0] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],
522 trs,primaryVertex[2],pflag);
523 nClusLay[1] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],
524 trs,primaryVertex[2],pflag);
527 fPoint3[0]=fPointc[0];
528 fPoint3[1]=fPointc[1];
530 nClusLay[2]= SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],
531 trs,primaryVertex[2],pflag);
536 nClusLay[3]= SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],
537 trs,primaryVertex[2],pflag);
542 nClusLay[4]=SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],
543 trs,primaryVertex[2],pflag);
546 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-1;nnp++){
547 if(nClusLay[nnp]!=0) layOK+=1;
549 if(layOK>=minNPoints){ // 5/6
550 //printf("--NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
552 AliITStrackV2* tr2 = 0;
553 tr2 = FitTrack(trs,primaryVertex);
555 //printf("--NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
557 StoreTrack(tr2,event);
563 }//end loop on clusters of layer2
568 // search for tracks starting from SPD2, SDD1, SDD2, SSD2
569 // for cosmics (A. Dainese 31.07.07)
570 if(fOuterStartLayer>0 && !AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad()) {
571 for(Int_t innLay=1; innLay<=fOuterStartLayer; innLay++) {
572 //printf("Searching from layer %d outward\n",innLay);
573 minNPoints=AliITSgeomTGeo::GetNLayers()-innLay;
574 for(Int_t i=innLay;i<AliITSgeomTGeo::GetNLayers();i++)
575 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i))
577 // counter for clusters on each layer
579 for(Int_t nloop=0;nloop<fNloop;nloop++){
580 Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
581 while(nclInnLay--){ //loop starting from layer innLay
584 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
587 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
588 fPhic = arr->GetPhi();
589 fLambdac = arr->GetLambda();
590 fPhiEstimate = fPhic;
592 AliITStrackSA* trs = new AliITStrackSA();
593 fPoint1[0]=primaryVertex[0];
594 fPoint1[1]=primaryVertex[1];
595 fPoint2[0]=arr->GetX();
596 fPoint2[1]=arr->GetY();
599 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
602 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
603 trs,primaryVertex[2],pflag);
604 for(Int_t nextLay=innLay+1; nextLay<AliITSgeomTGeo::GetNLayers(); nextLay++) {
606 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
607 trs,primaryVertex[2],pflag);
611 fPoint3[0]=fPointc[0];
612 fPoint3[1]=fPointc[1];
620 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
621 if(nClusLay[nnp]!=0) layOK+=1;
623 if(layOK>=minNPoints){
624 //printf("---NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
625 AliITStrackV2* tr2 = 0;
626 tr2 = FitTrack(trs,primaryVertex);
628 //printf("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
630 StoreTrack(tr2,event);
636 }//end loop on clusters of innLay
637 } //end loop on window sizes
639 } //end loop on innLay
640 } //end if(fOuterStartLayer>0)
643 // search for 1-point tracks in SPD, only for cosmics
644 // (A.Dainese 21.03.08)
645 if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
646 TMath::Abs(event->GetMagneticField())<0.01) {
647 Int_t outerLayer=1; // only SPD
648 for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
649 // counter for clusters on each layer
651 for(Int_t nloop=0;nloop<fNloop;nloop++){
652 Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
653 while(nclInnLay--){ //loop starting from layer innLay
656 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
659 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
660 fPhic = arr->GetPhi();
661 fLambdac = arr->GetLambda();
662 fPhiEstimate = fPhic;
664 AliITStrackSA* trs = new AliITStrackSA();
665 fPoint1[0]=primaryVertex[0];
666 fPoint1[1]=primaryVertex[1];
667 fPoint2[0]=arr->GetX();
668 fPoint2[1]=arr->GetY();
671 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
674 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
675 trs,primaryVertex[2],pflag);
676 for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
678 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
679 trs,primaryVertex[2],pflag);
683 fPoint3[0]=fPointc[0];
684 fPoint3[1]=fPointc[1];
692 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
693 if(nClusLay[nnp]!=0) layOK+=1;
696 //printf("----NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]);
697 AliITStrackV2* tr2 = 0;
698 Bool_t onePoint = kTRUE;
699 tr2 = FitTrack(trs,primaryVertex,onePoint);
701 //printf("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters());
703 StoreTrack(tr2,event);
709 }//end loop on clusters of innLay
710 } //end loop on window sizes
712 } //end loop on innLay
713 } // end search 1-point tracks
715 Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks());
720 //________________________________________________________________________
722 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
723 //fit of the found track (most general case, also <6 points, layers missing)
724 // A.Dainese 16.11.07
727 const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
729 static Int_t firstmod[AliITSgeomTGeo::kNLayers];
731 static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
732 static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
733 static Int_t nnn[AliITSgeomTGeo::kNLayers];
734 static Int_t kkk[AliITSgeomTGeo::kNLayers];
735 static Int_t end[AliITSgeomTGeo::kNLayers];
736 static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
738 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
739 firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
743 for(Int_t j=0;j<kMaxClu; j++){
751 Int_t nclusters = tr->GetNumberOfClustersSA();
752 for(Int_t ncl=0;ncl<nclusters;ncl++){
753 Int_t index = tr->GetClusterIndexSA(ncl);
754 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
755 if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
756 Int_t lay = (index & 0xf0000000) >> 28;
757 Int_t nInLay=end[lay];
758 listlayer[lay][nInLay]=cl;
760 clind[lay][ind]=index;
764 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
765 for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
766 Int_t mark = tr->GetClusterMark(nlay,ncl);
768 clmark[nlay][ind]=mark;
773 Int_t firstLay=-1,secondLay=-1;
774 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
780 } else if(secondLay==-1) {
786 if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
788 for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
789 AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0][l0];
790 for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
791 AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1][l1];
792 for(Int_t l2=0;l2<end[2];l2++){ //loop on layer 3
793 AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2][l2];
794 for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
795 AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3][l3];
796 for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
797 AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4][l4];
798 for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
799 AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5][l5];
802 Double_t x1,y1,z1,sx1,sy1,sz1;
803 Double_t x2,y2,z2,sx2,sy2,sz2;
804 AliITSRecPoint* p1=0;
805 AliITSRecPoint* p2=0;
806 Int_t index1=0,index2=0;
812 index1=clind[0][l0];mrk1=clmark[0][l0];
816 index1=clind[1][l1];mrk1=clmark[1][l1];
820 index1=clind[2][l2];mrk1=clmark[2][l2];
824 index1=clind[3][l3];mrk1=clmark[3][l3];
828 index1=clind[4][l4];mrk1=clmark[4][l4];
835 index2=clind[1][l1];mrk2=clmark[1][l1];
839 index2=clind[2][l2];mrk2=clmark[2][l2];
843 index2=clind[3][l3];mrk2=clmark[3][l3];
847 index2=clind[4][l4];mrk2=clmark[4][l4];
851 index2=clind[5][l5];mrk2=clmark[5][l5];
859 Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
860 Int_t layer,ladder,detector;
861 AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
862 Float_t yclu1 = p1->GetY();
863 Float_t zclu1 = p1->GetZ();
864 Double_t cv=0,tgl2=0,phi2=0;
867 AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
877 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,cln2);
884 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
885 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
886 phi2 = TMath::ATan2((y2-y1),(x2-x1));
887 } else { // special case of 1-point tracks, only for cosmics (B=0)
888 x2 = primaryVertex[0];
889 y2 = primaryVertex[1];
890 z2 = primaryVertex[2];
892 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
893 phi2 = TMath::ATan2((y1-y2),(x1-x2));
897 AliITStrackSA* trac = new AliITStrackSA(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
901 trac->AddClusterV2(5,(clind[5][l5] & 0x0fffffff)>>0);
902 trac->AddClusterMark(5,clmark[5][l5]);
905 trac->AddClusterV2(4,(clind[4][l4] & 0x0fffffff)>>0);
906 trac->AddClusterMark(4,clmark[4][l4]);
909 trac->AddClusterV2(3,(clind[3][l3] & 0x0fffffff)>>0);
910 trac->AddClusterMark(3,clmark[3][l3]);
913 trac->AddClusterV2(2,(clind[2][l2] & 0x0fffffff)>>0);
914 trac->AddClusterMark(2,clmark[2][l2]);
917 trac->AddClusterV2(1,(clind[1][l1] & 0x0fffffff)>>0);
918 trac->AddClusterMark(1,clmark[1][l1]);
921 trac->AddClusterV2(0,(clind[0][l0] & 0x0fffffff)>>0);
922 trac->AddClusterMark(0,clmark[0][l0]);
925 //fit with Kalman filter using AliITStrackerMI::RefitAt()
926 AliITStrackMI* ot = new AliITStrackSA(*trac);
928 ot->ResetCovariance(10.);
931 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),ot,trac)){ //fit from layer 1 to layer 6
932 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
933 otrack2->ResetCovariance(10.);
934 otrack2->ResetClusters();
935 //fit from layer 6 to layer 1
936 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),otrack2,ot)) {
937 fListOfTracks->AddLast(otrack2);
938 fListOfSATracks->AddLast(trac);
957 if(fListOfTracks->GetEntries()==0) return 0;
959 Int_t lowchi2 = FindTrackLowChiSquare();
960 AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
961 AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
963 if(otrack==0) return 0;
965 Int_t indexc[AliITSgeomTGeo::kNLayers];
966 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) indexc[i]=0;
967 for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
968 indexc[nind] = otrack->GetClusterIndex(nind);
971 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
972 if(i<otrack->GetNumberOfClusters()) {
973 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc[i]);
974 labl[i][0]=cl->GetLabel(0);
975 labl[i][1]=cl->GetLabel(1);
976 labl[i][2]=cl->GetLabel(2);
984 CookLabel(otrack,0.); //MI change - to see fake ratio
986 Int_t label=FindLabel(labl[0][0],labl[1][0],labl[2][0],labl[3][0],labl[4][0],labl[5][0]);
988 for(Int_t i=0;i<otrack->GetNumberOfClusters();i++)
989 if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
991 if(lflag<otrack->GetNumberOfClusters()) label = -label;
992 otrack->SetLabel(label);
994 //remove clusters of found track
995 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
996 for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
997 Int_t index = trsa->GetClusterMark(nlay,cln);
998 fCluLayer[nlay]->RemoveAt(index);
999 RemoveClusterCoord(nlay,index);
1000 fCluLayer[nlay]->Compress();
1008 //_______________________________________________________
1009 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event) const
1012 // Add new track to the ESD
1014 AliESDtrack outtrack;
1015 outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
1016 for(Int_t i=0;i<12;i++) {
1017 outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
1019 event->AddTrack(&outtrack);
1025 //_______________________________________________________
1026 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
1027 //function used to to find the clusters associated to the track
1029 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(layer)) return 0;
1032 AliITSlayer &lay = fgLayers[layer];
1033 Double_t r=lay.GetR();
1035 Float_t cx1,cx2,cy1,cy2;
1036 FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
1037 if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
1039 Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
1040 Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
1041 fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
1045 Int_t ncl = fCluLayer[layer]->GetEntries();
1046 for (Int_t index=0; index<ncl; index++) {
1047 AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
1049 if (c->GetQ()<=0) continue;
1050 if(layer>1 && c->GetQ()<=fMinQ) continue;
1052 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
1053 Double_t phi = arr->GetPhi();
1054 if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
1056 Double_t lambda = arr->GetLambda();
1057 if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
1059 if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
1060 if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
1061 Int_t orind = arr->GetOrInd();
1062 trs->AddClusterSA(layer,orind);
1063 trs->AddClusterMark(layer,index);
1069 fPointc[0]=arr->GetX();
1070 fPointc[1]=arr->GetY();
1076 //________________________________________________________________
1077 void AliITStrackerSA::UpdatePoints(){
1078 //update of points for the estimation of the curvature
1080 fPoint2[0]=fPoint3[0];
1081 fPoint2[1]=fPoint3[1];
1082 fPoint3[0]=fPointc[0];
1083 fPoint3[1]=fPointc[1];
1088 //___________________________________________________________________
1089 Int_t AliITStrackerSA::FindEquation(Float_t x1, Float_t y1, Float_t x2, Float_t y2, Float_t x3, Float_t y3,Float_t& a, Float_t& b, Float_t& c){
1091 //given (x,y) of three recpoints (in global coordinates)
1092 //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
1094 Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1095 if(den==0) return 0;
1096 a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1097 b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1098 c = -x1*x1-y1*y1-a*x1-b*y1;
1101 //__________________________________________________________________________
1102 Int_t AliITStrackerSA::FindIntersection(Float_t a1, Float_t b1, Float_t c1, Float_t c2,Float_t& x1,Float_t& y1, Float_t& x2, Float_t& y2){
1104 //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
1105 //c2 is -rlayer*rlayer
1109 Double_t aA = (b1*b1)/(a1*a1)+1;
1110 Double_t bB = (-2*m*b1/(a1*a1));
1111 Double_t cC = c2+(m*m)/(a1*a1);
1112 Double_t dD = bB*bB-4*aA*cC;
1115 y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
1116 y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
1117 x1 = (c2-c1-b1*y1)/a1;
1118 x2 = (c2-c1-b1*y2)/a1;
1122 //____________________________________________________________________
1123 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
1124 x2,Double_t y2,Double_t x3,Double_t y3){
1126 //calculates the curvature of track
1127 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1128 if(den==0) return 0;
1129 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1130 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1131 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
1134 if((a*a+b*b-4*c)<0) return 0;
1135 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
1136 if(rad==0) return 0;
1138 if((x1>0 && y1>0 && x1<xc)) rad*=-1;
1139 if((x1<0 && y1>0 && x1<xc)) rad*=-1;
1140 // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
1141 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
1148 //____________________________________________________________________
1149 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
1151 //Returns the point closest to pp
1153 Double_t diff1 = p1-pp;
1154 Double_t diff2 = p2-pp;
1156 if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
1157 else fPhiEstimate=p2;
1158 return fPhiEstimate;
1163 //_________________________________________________________________
1164 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
1165 // returns track with lowest chi square
1166 Int_t dim=fListOfTracks->GetEntries();
1167 if(dim<=1) return 0;
1168 AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
1169 Double_t minChi2=trk->GetChi2();
1171 for(Int_t i=1;i<dim;i++){
1172 trk = (AliITStrackV2*)fListOfTracks->At(i);
1173 Double_t chi2=trk->GetChi2();
1182 //__________________________________________________________
1183 Int_t AliITStrackerSA::FindLabel(Int_t l0, Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5){
1185 //function used to determine the track label
1187 Int_t lb[6] = {l0,l1,l2,l3,l4,l5};
1188 Int_t aa[6]={1,1,1,1,1,1};
1191 Int_t k=0;Int_t w=0;Int_t num=6;
1192 for(Int_t i=5;i>=0;i--) if(lb[i]==-1) num=i;
1196 for(Int_t i=k+1;i<num;i++){
1198 if(lb[k]==lb[i] && aa[k]!=0){
1209 for(Int_t j=0;j<6;j++){
1222 if(num<1) return -1;
1226 //_____________________________________________________________________________
1227 Int_t AliITStrackerSA::Label(Int_t gl1, Int_t gl2, Int_t gl3, Int_t gl4, Int_t gl5, Int_t gl6,Int_t gl7, Int_t gl8, Int_t gl9, Int_t gl10,Int_t gl11,
1228 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t minNPoints){
1231 //function used to assign label to the found track. If track is fake, the label is negative
1233 Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
1234 Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
1235 Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
1236 Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
1237 Int_t lflag=0;Int_t num=6;
1238 if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
1240 for(Int_t i=0;i<num;i++){
1241 if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
1244 if(lflag>=minNPoints) return ll;
1249 //_____________________________________________________________________________
1250 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
1251 // Set sizes of the phi and lambda windows used for track finding
1253 if(fPhiWin) delete [] fPhiWin;
1254 if(fLambdaWin) delete [] fLambdaWin;
1255 fPhiWin = new Double_t[fNloop];
1256 fLambdaWin = new Double_t[fNloop];
1257 Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
1258 Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
1259 for(Int_t k=0;k<fNloop;k++){
1260 Float_t phi=phimin+k*stepPhi;
1261 Float_t lam=lambdamin+k*stepLambda;
1266 //_____________________________________________________________________________
1267 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1268 // Set sizes of the phi and lambda windows used for track finding
1270 if(phi){ // user defined values
1271 fPhiWin = new Double_t[fNloop];
1272 fLambdaWin = new Double_t[fNloop];
1273 for(Int_t k=0;k<fNloop;k++){
1275 fLambdaWin[k]=lam[k];
1278 else { // default values
1280 Double_t phid[33] = {0.002,0.003,0.004,0.0045,0.0047,
1281 0.005,0.0053,0.0055,
1282 0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1283 0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1284 0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1285 Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1287 0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1288 0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1289 0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1296 fPhiWin = new Double_t[fNloop];
1297 fLambdaWin = new Double_t[fNloop];
1299 Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1301 for(Int_t k=0;k<fNloop;k++){
1302 fPhiWin[k]=phid[k]*factor;
1303 fLambdaWin[k]=lambdad[k]*factor;
1309 //_______________________________________________________________________
1310 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1311 //Returns values of phi (azimuthal) and lambda angles for a given cluster
1313 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1314 Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1315 Float_t tx,ty,tz; fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1317 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1318 Double_t phi1=TMath::Pi()/2+alpha;
1319 if (lay==1) phi1+=TMath::Pi();
1321 Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1322 Float_t r=tx*cp+ty*sp;
1324 xyz= r*cp - cl->GetY()*sp;
1325 y= r*sp + cl->GetY()*cp;
1329 cl->GetGlobalXYZ(xyz);
1334 phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1335 lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1338 //________________________________________________________________________
1339 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1341 //returns sigmax, y, z of cluster in global coordinates
1343 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1345 AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1347 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1348 Double_t phi=TMath::Pi()/2+alpha;
1349 if (lay==1) phi+=TMath::Pi();
1351 Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1354 cl->GetGlobalCov(covm);
1355 sx=TMath::Sqrt(covm[0]);
1356 sy=TMath::Sqrt(covm[3]);
1357 sz=TMath::Sqrt(covm[5]);
1359 sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1360 sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1361 sz = TMath::Sqrt(cl->GetSigmaZ2());