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"
45 ClassImp(AliITStrackerSA)
47 //____________________________________________________________________________
48 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
71 // Default constructor
75 //____________________________________________________________________________
76 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
100 // Standard constructor (Vertex is known and passed to this obj.)
102 AliWarning("\"geom\" is actually a dummy argument !");
110 //____________________________________________________________________________
111 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliESDVertex *vert):AliITStrackerMI(0),
135 // Standard constructor (Vertex is known and passed to this obj.)
137 AliWarning("\"geom\" is actually a dummy argument !");
143 //____________________________________________________________________________
144 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliITSVertexer *vertexer):AliITStrackerMI(0),
168 // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
170 AliWarning("\"geom\" is actually a dummy argument !");
173 fVertexer = vertexer;
177 //____________________________________________________________________________
178 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(),
179 fPhiEstimate(tracker.fPhiEstimate),
180 fITSStandAlone(tracker.fITSStandAlone),
181 fLambdac(tracker.fLambdac),
182 fPhic(tracker.fPhic),
183 fCoef1(tracker.fCoef1),
184 fCoef2(tracker.fCoef2),
185 fCoef3(tracker.fCoef3),
186 fNloop(tracker.fNloop),
187 fPhiWin(tracker.fPhiWin),
188 fLambdaWin(tracker.fLambdaWin),
189 fVert(tracker.fVert),
190 fVertexer(tracker.fVertexer),
191 fListOfTracks(tracker.fListOfTracks),
192 fListOfSATracks(tracker.fListOfSATracks),
193 fITSclusters(tracker.fITSclusters),
194 fInwardFlag(tracker.fInwardFlag),
195 fOuterStartLayer(tracker.fOuterStartLayer),
196 fInnerStartLayer(tracker.fInnerStartLayer),
197 fMinNPoints(tracker.fMinNPoints),
198 fMinQ(tracker.fMinQ),
199 fCluLayer(tracker.fCluLayer),
200 fCluCoord(tracker.fCluCoord) {
202 for(Int_t i=0;i<2;i++){
203 fPoint1[i]=tracker.fPoint1[i];
204 fPoint2[i]=tracker.fPoint2[i];
205 fPoint3[i]=tracker.fPoint3[i];
206 fPointc[i]=tracker.fPointc[i];
208 if(tracker.fVertexer && tracker.fVert){
209 fVert = new AliESDVertex(*tracker.fVert);
212 fVert = tracker.fVert;
214 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
215 fCluLayer[i] = tracker.fCluLayer[i];
216 fCluCoord[i] = tracker.fCluCoord[i];
219 //______________________________________________________________________
220 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& source){
221 // Assignment operator.
222 this->~AliITStrackerSA();
223 new(this) AliITStrackerSA(source);
228 //____________________________________________________________________________
229 AliITStrackerSA::~AliITStrackerSA(){
231 // if fVertexer is not null, the AliESDVertex obj. is owned by this class
232 // and is deleted here
234 if(fVert)delete fVert;
239 if(fPhiWin)delete []fPhiWin;
240 if(fLambdaWin)delete []fLambdaWin;
241 fListOfTracks->Delete();
242 delete fListOfTracks;
243 fListOfSATracks->Delete();
244 delete fListOfSATracks;
246 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
248 fCluLayer[i]->Delete();
255 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
257 fCluCoord[i]->Delete();
266 //____________________________________________________________________________
267 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
268 // This method is used to find and fit the tracks. By default the corresponding
269 // method in the parent class is invoked. In this way a combined tracking
270 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
271 // is done in the ITS only. In the standard reconstruction chain this option
272 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
275 rc=AliITStrackerMI::Clusters2Tracks(event);
278 AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
280 if(!rc) rc=FindTracks(event);
284 //____________________________________________________________________________
285 void AliITStrackerSA::Init(){
286 // Reset all data members
288 for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
298 Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
300 SetFixedWindowSizes();
302 Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
303 Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
304 Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
305 Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
306 SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
308 fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
310 SetOuterStartLayer(1);
312 fListOfTracks=new TObjArray(0,0);
313 fListOfSATracks=new TObjArray(0,0);
318 //_______________________________________________________________________
319 void AliITStrackerSA::ResetForFinding(){
320 // Reset data members used in all loops during track finding
322 for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
330 fListOfTracks->Delete();
331 fListOfSATracks->Delete();
336 //______________________________________________________________________
337 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event){
339 // Track finder using the ESD object
341 AliDebug(2,Form(" field is %f",event->GetMagneticField()));
342 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
345 Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
348 //Reads event and mark clusters of traks already found, with flag kITSin
349 Int_t nentr=event->GetNumberOfTracks();
350 if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kFALSE) {
352 AliESDtrack *track=event->GetTrack(nentr);
353 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
355 Int_t ncl = track->GetITSclusters(idx);
356 for(Int_t k=0;k<ncl;k++){
357 AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
358 cll->SetBit(kSAflag);
364 Double_t primaryVertex[3];
365 event->GetVertex()->GetXYZ(primaryVertex);
366 //Creates TClonesArray with clusters for each layer. The clusters already used
367 //by AliITStrackerMI are not considered
368 Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
369 Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
370 if (fCluLayer == 0) {
371 fCluLayer = new TClonesArray*[AliITSgeomTGeo::kNLayers];
372 fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
373 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
378 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
379 AliITSlayer &layer=fgLayers[i];
380 if (!ForceSkippingOfLayer(i)) {
381 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
382 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
383 if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
384 if(cls->GetQ()==0) continue; //fake clusters dead zones
385 if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
391 fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
393 fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
396 for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
397 TClonesArray &clulay = *fCluLayer[ilay];
398 TClonesArray &clucoo = *fCluCoord[ilay];
399 AliITSlayer &layer=fgLayers[ilay];
400 if (!ForceSkippingOfLayer(ilay)) {
401 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
402 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
403 if(cls->TestBit(kSAflag)==kTRUE) continue;
404 if(cls->GetQ()==0) continue;
405 if(ilay>1 && cls->GetQ()<=fMinQ) continue;
406 Double_t phi=0;Double_t lambda=0;
407 Float_t x=0;Float_t y=0;Float_t z=0;
408 Float_t sx=0;Float_t sy=0;Float_t sz=0;
409 GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
410 GetCoorErrors(cls,sx,sy,sz);
411 new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
412 new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
421 static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
423 // loop on minimum number of points
424 for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) {
426 if(!fInwardFlag){ // Tracking outwards from the inner layers
427 // loop on starting layer for track finding
428 for(Int_t innLay=0; innLay<=fOuterStartLayer; innLay++) {
429 if(ForceSkippingOfLayer(innLay)) continue;
430 Int_t minNPoints=iMinNPoints-innLay;
431 for(Int_t i=innLay+1;i<AliITSgeomTGeo::GetNLayers();i++)
432 if(ForceSkippingOfLayer(i))
434 if(minNPoints<fMinNPoints) continue;
436 // loop on phi and lambda window size
437 for(Int_t nloop=0;nloop<fNloop;nloop++){
438 Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
441 Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
443 AliITStrackSA* trs = new AliITStrackSA();
447 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
450 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
451 trs,primaryVertex[2],pflag);
452 for(Int_t nextLay=innLay+1; nextLay<AliITSgeomTGeo::GetNLayers(); nextLay++) {
454 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
455 trs,primaryVertex[2],pflag);
459 fPoint3[0]=fPointc[0];
460 fPoint3[1]=fPointc[1];
468 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
469 if(nClusLay[nnp]!=0) layOK+=1;
471 if(layOK>=minNPoints){
472 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]));
473 AliITStrackV2* tr2 = 0;
474 tr2 = FitTrack(trs,primaryVertex);
476 AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
478 StoreTrack(tr2,event);
484 }//end loop on clusters of innLay
485 } //end loop on window sizes
486 } //end loop on innLay
487 }else{// Tracking inwards from the outer layers
488 // loop on starting layer for track finding
489 for(Int_t outLay=AliITSgeomTGeo::GetNLayers()-1; outLay>=fInnerStartLayer; outLay--) {
491 if(ForceSkippingOfLayer(outLay)) continue;
492 Int_t minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-outLay);
493 for(Int_t i=0;i<outLay-1;i++)
494 if(ForceSkippingOfLayer(i))
496 if(minNPoints<fMinNPoints) continue;
498 // loop on phi and lambda window size
499 for(Int_t nloop=0;nloop<fNloop;nloop++){
500 Int_t nclOutLay=fCluLayer[outLay]->GetEntries();
503 Bool_t useRP=SetFirstPoint(outLay,nclOutLay,primaryVertex);
505 AliITStrackSA* trs = new AliITStrackSA();
509 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
512 nClusLay[kk] = SearchClusters(outLay,fPhiWin[nloop],fLambdaWin[nloop],
513 trs,primaryVertex[2],pflag);
514 for(Int_t nextLay=outLay-1; nextLay>=0; nextLay--) {
516 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
517 trs,primaryVertex[2],pflag);
521 fPoint3[0]=fPointc[0];
522 fPoint3[1]=fPointc[1];
530 for(Int_t nnp=outLay; nnp>=0; nnp--){
531 if(nClusLay[nnp]!=0) layOK+=1;
533 if(layOK>=minNPoints){
534 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]));
535 AliITStrackV2* tr2 = 0;
536 tr2 = FitTrack(trs,primaryVertex);
538 AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
540 StoreTrack(tr2,event);
546 }//end loop on clusters of innLay
547 } //end loop on window sizes
548 } //end loop on innLay
549 } // end if (fInwardFlag)
550 }//end loop on min points
552 // search for 1-point tracks in SPD, only for cosmics
553 // (A.Dainese 21.03.08)
554 if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
555 TMath::Abs(event->GetMagneticField())<0.01) {
556 Int_t outerLayer=1; // only SPD
557 for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
558 // counter for clusters on each layer
560 for(Int_t nloop=0;nloop<fNloop;nloop++){
561 Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
562 while(nclInnLay--){ //loop starting from layer innLay
564 Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
566 AliITStrackSA* trs = new AliITStrackSA();
570 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
573 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
574 trs,primaryVertex[2],pflag);
575 for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
577 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
578 trs,primaryVertex[2],pflag);
582 fPoint3[0]=fPointc[0];
583 fPoint3[1]=fPointc[1];
591 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
592 if(nClusLay[nnp]!=0) layOK+=1;
595 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]));
596 AliITStrackV2* tr2 = 0;
597 Bool_t onePoint = kTRUE;
598 tr2 = FitTrack(trs,primaryVertex,onePoint);
600 AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
602 StoreTrack(tr2,event);
608 }//end loop on clusters of innLay
609 } //end loop on window sizes
611 } //end loop on innLay
612 } // end search 1-point tracks
614 AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
620 //________________________________________________________________________
622 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
623 //fit of the found track (most general case, also <6 points, layers missing)
624 // A.Dainese 16.11.07
627 const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
629 static Int_t firstmod[AliITSgeomTGeo::kNLayers];
631 static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
632 static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
633 static Int_t end[AliITSgeomTGeo::kNLayers];
634 static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
636 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
637 firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
639 for(Int_t j=0;j<kMaxClu; j++){
647 Int_t nclusters = tr->GetNumberOfClustersSA();
648 for(Int_t ncl=0;ncl<nclusters;ncl++){
649 Int_t index = tr->GetClusterIndexSA(ncl);
650 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
651 if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
652 Int_t lay = (index & 0xf0000000) >> 28;
653 Int_t nInLay=end[lay];
654 listlayer[lay][nInLay]=cl;
655 clind[lay][nInLay]=index;
659 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
660 for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
661 Int_t mark = tr->GetClusterMark(nlay,ncl);
662 clmark[nlay][ncl]=mark;
667 Int_t firstLay=-1,secondLay=-1;
668 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
674 } else if(secondLay==-1) {
680 if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
682 for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
683 AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0][l0];
684 for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
685 AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1][l1];
686 for(Int_t l2=0;l2<end[2];l2++){ //loop on layer 3
687 AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2][l2];
688 for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
689 AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3][l3];
690 for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
691 AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4][l4];
692 for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
693 AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5][l5];
696 Double_t x1,y1,z1,sx1,sy1,sz1;
697 Double_t x2,y2,z2,sx2,sy2,sz2;
698 AliITSRecPoint* p1=0;
699 AliITSRecPoint* p2=0;
700 Int_t index1=0,index2=0;
706 index1=clind[0][l0];mrk1=clmark[0][l0];
710 index1=clind[1][l1];mrk1=clmark[1][l1];
714 index1=clind[2][l2];mrk1=clmark[2][l2];
718 index1=clind[3][l3];mrk1=clmark[3][l3];
722 index1=clind[4][l4];mrk1=clmark[4][l4];
729 index2=clind[1][l1];mrk2=clmark[1][l1];
733 index2=clind[2][l2];mrk2=clmark[2][l2];
737 index2=clind[3][l3];mrk2=clmark[3][l3];
741 index2=clind[4][l4];mrk2=clmark[4][l4];
745 index2=clind[5][l5];mrk2=clmark[5][l5];
753 Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
754 Int_t layer,ladder,detector;
755 AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
756 Float_t yclu1 = p1->GetY();
757 Float_t zclu1 = p1->GetZ();
758 Double_t cv=0,tgl2=0,phi2=0;
761 AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
771 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,cln2);
778 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
779 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
780 phi2 = TMath::ATan2((y2-y1),(x2-x1));
781 } else { // special case of 1-point tracks, only for cosmics (B=0)
782 x2 = primaryVertex[0];
783 y2 = primaryVertex[1];
784 z2 = primaryVertex[2];
786 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
787 phi2 = TMath::ATan2((y1-y2),(x1-x2));
791 AliITStrackSA* trac = new AliITStrackSA(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
795 trac->AddClusterV2(5,(clind[5][l5] & 0x0fffffff)>>0);
796 trac->AddClusterMark(5,clmark[5][l5]);
799 trac->AddClusterV2(4,(clind[4][l4] & 0x0fffffff)>>0);
800 trac->AddClusterMark(4,clmark[4][l4]);
803 trac->AddClusterV2(3,(clind[3][l3] & 0x0fffffff)>>0);
804 trac->AddClusterMark(3,clmark[3][l3]);
807 trac->AddClusterV2(2,(clind[2][l2] & 0x0fffffff)>>0);
808 trac->AddClusterMark(2,clmark[2][l2]);
811 trac->AddClusterV2(1,(clind[1][l1] & 0x0fffffff)>>0);
812 trac->AddClusterMark(1,clmark[1][l1]);
815 trac->AddClusterV2(0,(clind[0][l0] & 0x0fffffff)>>0);
816 trac->AddClusterMark(0,clmark[0][l0]);
819 //fit with Kalman filter using AliITStrackerMI::RefitAt()
820 AliITStrackMI* ot = new AliITStrackSA(*trac);
822 ot->ResetCovariance(10.);
825 // Propagate inside the innermost layer with a cluster
826 if(ot->Propagate(ot->GetX()-0.1*ot->GetX())) {
828 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),ot,trac)){ //fit from layer 1 to layer 6
829 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
830 otrack2->ResetCovariance(10.);
831 otrack2->ResetClusters();
832 //fit from layer 6 to layer 1
833 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),otrack2,ot)) {
834 fListOfTracks->AddLast(otrack2);
835 fListOfSATracks->AddLast(trac);
854 if(fListOfTracks->GetEntries()==0) return 0;
856 Int_t lowchi2 = FindTrackLowChiSquare();
857 AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
858 AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
860 if(otrack==0) return 0;
862 Int_t indexc[AliITSgeomTGeo::kNLayers];
863 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) indexc[i]=0;
864 for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
865 indexc[nind] = otrack->GetClusterIndex(nind);
868 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
869 if(i<otrack->GetNumberOfClusters()) {
870 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc[i]);
871 labl[i][0]=cl->GetLabel(0);
872 labl[i][1]=cl->GetLabel(1);
873 labl[i][2]=cl->GetLabel(2);
881 CookLabel(otrack,0.); //MI change - to see fake ratio
883 Int_t label=FindLabel(labl[0][0],labl[1][0],labl[2][0],labl[3][0],labl[4][0],labl[5][0]);
885 for(Int_t i=0;i<otrack->GetNumberOfClusters();i++)
886 if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
888 if(lflag<otrack->GetNumberOfClusters()) label = -label;
889 otrack->SetLabel(label);
891 //remove clusters of found track
892 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
893 for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
894 Int_t index = trsa->GetClusterMark(nlay,cln);
895 fCluLayer[nlay]->RemoveAt(index);
896 RemoveClusterCoord(nlay,index);
897 fCluLayer[nlay]->Compress();
905 //_______________________________________________________
906 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event) const
909 // Add new track to the ESD
911 AliESDtrack outtrack;
912 outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
913 for(Int_t i=0;i<12;i++) {
914 outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
916 Double_t sdedx[4]={0.,0.,0.,0.};
917 for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
918 outtrack.SetITSdEdxSamples(sdedx);
919 event->AddTrack(&outtrack);
925 //_______________________________________________________
926 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
927 //function used to to find the clusters associated to the track
929 if(ForceSkippingOfLayer(layer)) return 0;
932 AliITSlayer &lay = fgLayers[layer];
933 Double_t r=lay.GetR();
935 Float_t cx1,cx2,cy1,cy2;
936 FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
937 if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
939 Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
940 Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
941 fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
945 Int_t ncl = fCluLayer[layer]->GetEntries();
946 for (Int_t index=0; index<ncl; index++) {
947 AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
949 if (c->GetQ()<=0) continue;
950 if(layer>1 && c->GetQ()<=fMinQ) continue;
952 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
953 Double_t phi = arr->GetPhi();
954 if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
956 Double_t lambda = arr->GetLambda();
957 if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
959 if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
960 if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
961 Int_t orind = arr->GetOrInd();
962 trs->AddClusterSA(layer,orind);
963 trs->AddClusterMark(layer,index);
969 fPointc[0]=arr->GetX();
970 fPointc[1]=arr->GetY();
976 //________________________________________________________________
977 Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
978 // Sets the first point (seed) for tracking
980 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[lay]->At(clu);
981 if(!cl) return kFALSE;
982 if (cl->GetQ()<=0) return kFALSE;
983 if(lay>1 && cl->GetQ()<=fMinQ) return kFALSE;
985 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
986 fPhic = arr->GetPhi();
987 fLambdac = arr->GetLambda();
988 fPhiEstimate = fPhic;
989 fPoint1[0]=primaryVertex[0];
990 fPoint1[1]=primaryVertex[1];
991 fPoint2[0]=arr->GetX();
992 fPoint2[1]=arr->GetY();
996 //________________________________________________________________
997 void AliITStrackerSA::UpdatePoints(){
998 //update of points for the estimation of the curvature
1000 fPoint2[0]=fPoint3[0];
1001 fPoint2[1]=fPoint3[1];
1002 fPoint3[0]=fPointc[0];
1003 fPoint3[1]=fPointc[1];
1008 //___________________________________________________________________
1009 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){
1011 //given (x,y) of three recpoints (in global coordinates)
1012 //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
1014 Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1015 if(den==0) return 0;
1016 a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1017 b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1018 c = -x1*x1-y1*y1-a*x1-b*y1;
1021 //__________________________________________________________________________
1022 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){
1024 //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
1025 //c2 is -rlayer*rlayer
1029 Double_t aA = (b1*b1)/(a1*a1)+1;
1030 Double_t bB = (-2*m*b1/(a1*a1));
1031 Double_t cC = c2+(m*m)/(a1*a1);
1032 Double_t dD = bB*bB-4*aA*cC;
1035 y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
1036 y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
1037 x1 = (c2-c1-b1*y1)/a1;
1038 x2 = (c2-c1-b1*y2)/a1;
1042 //____________________________________________________________________
1043 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
1044 x2,Double_t y2,Double_t x3,Double_t y3){
1046 //calculates the curvature of track
1047 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1048 if(den==0) return 0;
1049 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1050 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1051 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
1054 if((a*a+b*b-4*c)<0) return 0;
1055 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
1056 if(rad==0) return 0;
1058 if((x1>0 && y1>0 && x1<xc)) rad*=-1;
1059 if((x1<0 && y1>0 && x1<xc)) rad*=-1;
1060 // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
1061 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
1068 //____________________________________________________________________
1069 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
1071 //Returns the point closest to pp
1073 Double_t diff1 = p1-pp;
1074 Double_t diff2 = p2-pp;
1076 if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
1077 else fPhiEstimate=p2;
1078 return fPhiEstimate;
1083 //_________________________________________________________________
1084 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
1085 // returns track with lowest chi square
1086 Int_t dim=fListOfTracks->GetEntries();
1087 if(dim<=1) return 0;
1088 AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
1089 Double_t minChi2=trk->GetChi2();
1091 for(Int_t i=1;i<dim;i++){
1092 trk = (AliITStrackV2*)fListOfTracks->At(i);
1093 Double_t chi2=trk->GetChi2();
1102 //__________________________________________________________
1103 Int_t AliITStrackerSA::FindLabel(Int_t l0, Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5){
1105 //function used to determine the track label
1107 Int_t lb[6] = {l0,l1,l2,l3,l4,l5};
1108 Int_t aa[6]={1,1,1,1,1,1};
1111 Int_t k=0;Int_t w=0;Int_t num=6;
1112 for(Int_t i=5;i>=0;i--) if(lb[i]==-1) num=i;
1116 for(Int_t i=k+1;i<num;i++){
1118 if(lb[k]==lb[i] && aa[k]!=0){
1129 for(Int_t j=0;j<6;j++){
1142 if(num<1) return -1;
1146 //_____________________________________________________________________________
1147 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,
1148 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t minNPoints){
1151 //function used to assign label to the found track. If track is fake, the label is negative
1153 Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
1154 Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
1155 Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
1156 Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
1157 Int_t lflag=0;Int_t num=6;
1158 if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
1160 for(Int_t i=0;i<num;i++){
1161 if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
1164 if(lflag>=minNPoints) return ll;
1169 //_____________________________________________________________________________
1170 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
1171 // Set sizes of the phi and lambda windows used for track finding
1173 if(fPhiWin) delete [] fPhiWin;
1174 if(fLambdaWin) delete [] fLambdaWin;
1175 fPhiWin = new Double_t[fNloop];
1176 fLambdaWin = new Double_t[fNloop];
1177 Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
1178 Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
1179 for(Int_t k=0;k<fNloop;k++){
1180 Float_t phi=phimin+k*stepPhi;
1181 Float_t lam=lambdamin+k*stepLambda;
1186 //_____________________________________________________________________________
1187 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1188 // Set sizes of the phi and lambda windows used for track finding
1190 if(phi){ // user defined values
1191 fPhiWin = new Double_t[fNloop];
1192 fLambdaWin = new Double_t[fNloop];
1193 for(Int_t k=0;k<fNloop;k++){
1195 fLambdaWin[k]=lam[k];
1198 else { // default values
1200 Double_t phid[33] = {0.002,0.003,0.004,0.0045,0.0047,
1201 0.005,0.0053,0.0055,
1202 0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1203 0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1204 0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1205 Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1207 0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1208 0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1209 0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1216 fPhiWin = new Double_t[fNloop];
1217 fLambdaWin = new Double_t[fNloop];
1219 Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1221 for(Int_t k=0;k<fNloop;k++){
1222 fPhiWin[k]=phid[k]*factor;
1223 fLambdaWin[k]=lambdad[k]*factor;
1229 //_______________________________________________________________________
1230 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1231 //Returns values of phi (azimuthal) and lambda angles for a given cluster
1233 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1234 Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1235 Float_t tx,ty,tz; fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1237 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1238 Double_t phi1=TMath::Pi()/2+alpha;
1239 if (lay==1) phi1+=TMath::Pi();
1241 Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1242 Float_t r=tx*cp+ty*sp;
1244 xyz= r*cp - cl->GetY()*sp;
1245 y= r*sp + cl->GetY()*cp;
1249 cl->GetGlobalXYZ(xyz);
1254 phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1255 lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1258 //________________________________________________________________________
1259 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1261 //returns sigmax, y, z of cluster in global coordinates
1263 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1265 AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1267 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1268 Double_t phi=TMath::Pi()/2+alpha;
1269 if (lay==1) phi+=TMath::Pi();
1271 Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1274 cl->GetGlobalCov(covm);
1275 sx=TMath::Sqrt(covm[0]);
1276 sy=TMath::Sqrt(covm[3]);
1277 sz=TMath::Sqrt(covm[5]);
1279 sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1280 sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1281 sz = TMath::Sqrt(cl->GetSigmaZ2());