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 ITS tracker class //
20 // Origin: Elisabetta Crescio - crescio@to.infn.it //
21 // Updated: Francesco Prino - prino@to.infn.it //
22 ///////////////////////////////////////////////////////////
28 #include <TObjArray.h>
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"
44 ClassImp(AliITStrackerSA)
46 //____________________________________________________________________________
47 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
70 // Default constructor
74 //____________________________________________________________________________
75 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
99 // Standard constructor (Vertex is known and passed to this obj.)
101 AliWarning("\"geom\" is actually a dummy argument !");
109 //____________________________________________________________________________
110 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliESDVertex *vert):AliITStrackerMI(0),
134 // Standard constructor (Vertex is known and passed to this obj.)
136 AliWarning("\"geom\" is actually a dummy argument !");
142 //____________________________________________________________________________
143 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliITSVertexer *vertexer):AliITStrackerMI(0),
167 // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
169 AliWarning("\"geom\" is actually a dummy argument !");
172 fVertexer = vertexer;
176 //____________________________________________________________________________
177 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(),
178 fPhiEstimate(tracker.fPhiEstimate),
179 fITSStandAlone(tracker.fITSStandAlone),
180 fLambdac(tracker.fLambdac),
181 fPhic(tracker.fPhic),
182 fCoef1(tracker.fCoef1),
183 fCoef2(tracker.fCoef2),
184 fCoef3(tracker.fCoef3),
185 fNloop(tracker.fNloop),
186 fPhiWin(tracker.fPhiWin),
187 fLambdaWin(tracker.fLambdaWin),
188 fVert(tracker.fVert),
189 fVertexer(tracker.fVertexer),
190 fListOfTracks(tracker.fListOfTracks),
191 fListOfSATracks(tracker.fListOfSATracks),
192 fITSclusters(tracker.fITSclusters),
193 fInwardFlag(tracker.fInwardFlag),
194 fOuterStartLayer(tracker.fOuterStartLayer),
195 fInnerStartLayer(tracker.fInnerStartLayer),
196 fMinNPoints(tracker.fMinNPoints),
197 fMinQ(tracker.fMinQ),
198 fCluLayer(tracker.fCluLayer),
199 fCluCoord(tracker.fCluCoord) {
201 for(Int_t i=0;i<2;i++){
202 fPoint1[i]=tracker.fPoint1[i];
203 fPoint2[i]=tracker.fPoint2[i];
204 fPoint3[i]=tracker.fPoint3[i];
205 fPointc[i]=tracker.fPointc[i];
207 if(tracker.fVertexer && tracker.fVert){
208 fVert = new AliESDVertex(*tracker.fVert);
211 fVert = tracker.fVert;
213 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
214 fCluLayer[i] = tracker.fCluLayer[i];
215 fCluCoord[i] = tracker.fCluCoord[i];
218 //______________________________________________________________________
219 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& source){
220 // Assignment operator.
221 this->~AliITStrackerSA();
222 new(this) AliITStrackerSA(source);
227 //____________________________________________________________________________
228 AliITStrackerSA::~AliITStrackerSA(){
230 // if fVertexer is not null, the AliESDVertex obj. is owned by this class
231 // and is deleted here
233 if(fVert)delete fVert;
238 if(fPhiWin)delete []fPhiWin;
239 if(fLambdaWin)delete []fLambdaWin;
240 fListOfTracks->Delete();
241 delete fListOfTracks;
242 fListOfSATracks->Delete();
243 delete fListOfSATracks;
245 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
247 fCluLayer[i]->Delete();
254 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
256 fCluCoord[i]->Delete();
265 //____________________________________________________________________________
266 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
267 // This method is used to find and fit the tracks. By default the corresponding
268 // method in the parent class is invoked. In this way a combined tracking
269 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
270 // is done in the ITS only. In the standard reconstruction chain this option
271 // 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");
281 rc=FindTracks(event,kFALSE);
283 const AliESDVertex *spdv = event->GetPrimaryVertexSPD();
284 if(spdv) nSPDcontr = spdv->GetNContributors();
285 if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE &&
286 nSPDcontr<=AliITSReconstructor::GetRecoParam()->GetMaxSPDcontrForSAToUseAllClusters()) {
287 rc=FindTracks(event,kTRUE);
293 //____________________________________________________________________________
294 void AliITStrackerSA::Init(){
295 // Reset all data members
297 for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
307 Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
309 SetFixedWindowSizes();
311 Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
312 Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
313 Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
314 Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
315 SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
317 fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
319 SetOuterStartLayer(1);
321 fListOfTracks=new TClonesArray("AliITStrackMI",100);
322 fListOfSATracks=new TClonesArray("AliITStrackSA",100);
327 //_______________________________________________________________________
328 void AliITStrackerSA::ResetForFinding(){
329 // Reset data members used in all loops during track finding
331 for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
339 fListOfTracks->Clear();
340 fListOfSATracks->Clear();
345 //______________________________________________________________________
346 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event, Bool_t useAllClusters){
348 // Track finder using the ESD object
350 AliDebug(2,Form(" field is %f",event->GetMagneticField()));
351 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
354 Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
357 //Reads event and mark clusters of traks already found, with flag kITSin
358 Int_t nentr=event->GetNumberOfTracks();
359 if(!useAllClusters) {
361 AliESDtrack *track=event->GetTrack(nentr);
362 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
364 Int_t ncl = track->GetITSclusters(idx);
365 for(Int_t k=0;k<ncl;k++){
366 AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
367 cll->SetBit(kSAflag);
373 AliESDtrack *track=event->GetTrack(nentr);
374 if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
376 Int_t ncl = track->GetITSclusters(idx);
377 for(Int_t k=0;k<ncl;k++){
378 AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
379 cll->ResetBit(kSAflag);
385 Double_t primaryVertex[3];
386 event->GetVertex()->GetXYZ(primaryVertex);
387 //Creates TClonesArray with clusters for each layer. The clusters already used
388 //by AliITStrackerMI are not considered
389 Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
390 Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
391 if (fCluLayer == 0) {
392 fCluLayer = new TClonesArray*[AliITSgeomTGeo::kNLayers];
393 fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
394 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
399 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
400 AliITSlayer &layer=fgLayers[i];
401 if (!ForceSkippingOfLayer(i)) {
402 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
403 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
404 if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
405 if(cls->GetQ()==0) continue; //fake clusters dead zones
406 if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
412 fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
414 fCluLayer[i]->Delete();
415 fCluLayer[i]->Expand(nclusters[i]);
418 fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
420 fCluCoord[i]->Delete();
421 fCluCoord[i]->Expand(nclusters[i]);
425 for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
426 TClonesArray &clulay = *fCluLayer[ilay];
427 TClonesArray &clucoo = *fCluCoord[ilay];
428 AliITSlayer &layer=fgLayers[ilay];
429 if (!ForceSkippingOfLayer(ilay)) {
430 for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
431 AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
432 if(cls->TestBit(kSAflag)==kTRUE) continue;
433 if(cls->GetQ()==0) continue;
434 if(ilay>1 && cls->GetQ()<=fMinQ) continue;
435 Double_t phi=0;Double_t lambda=0;
436 Float_t x=0;Float_t y=0;Float_t z=0;
437 Float_t sx=0;Float_t sy=0;Float_t sz=0;
438 GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
439 GetCoorErrors(cls,sx,sy,sz);
440 new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
441 new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
450 static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
451 Int_t startLayForSeed=0;
452 Int_t lastLayForSeed=fOuterStartLayer;
453 Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
456 startLayForSeed=AliITSgeomTGeo::GetNLayers()-1;
457 lastLayForSeed=fInnerStartLayer;
458 nSeedSteps=startLayForSeed-lastLayForSeed;
462 // loop on minimum number of points
463 for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) {
465 // loop on starting layer for track finding
466 for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
467 Int_t theLay=startLayForSeed+iSeedLay*seedStep;
468 if(ForceSkippingOfLayer(theLay)) continue;
469 Int_t minNPoints=iMinNPoints-theLay;
470 if(fInwardFlag) minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-theLay);
471 for(Int_t i=theLay+1;i<AliITSgeomTGeo::GetNLayers();i++)
472 if(ForceSkippingOfLayer(i))
474 if(minNPoints<fMinNPoints) continue;
476 // loop on phi and lambda window size
477 for(Int_t nloop=0;nloop<fNloop;nloop++){
478 Int_t nclTheLay=fCluLayer[theLay]->GetEntries();
481 Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
487 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
490 nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
491 &trs,primaryVertex[2],pflag);
492 Int_t nextLay=theLay+seedStep;
496 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
497 &trs,primaryVertex[2],pflag);
501 fPoint3[0]=fPointc[0];
502 fPoint3[1]=fPointc[1];
508 if(nextLay<0 || nextLay==6) goon=kFALSE;
514 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-theLay;nnp++){
515 if(nClusLay[nnp]!=0) layOK+=1;
518 for(Int_t nnp=theLay; nnp>=0; nnp--){
519 if(nClusLay[nnp]!=0) layOK+=1;
522 if(layOK>=minNPoints){
523 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]));
524 AliITStrackV2* tr2 = 0;
525 tr2 = FitTrack(&trs,primaryVertex);
529 AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
531 StoreTrack(tr2,event,useAllClusters);
536 }//end loop on clusters of theLay
537 } //end loop on window sizes
538 } //end loop on theLay
539 }//end loop on min points
541 // search for 1-point tracks in SPD, only for cosmics
542 // (A.Dainese 21.03.08)
543 if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
544 TMath::Abs(event->GetMagneticField())<0.01) {
545 Int_t outerLayer=1; // only SPD
546 for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
547 // counter for clusters on each layer
549 for(Int_t nloop=0;nloop<fNloop;nloop++){
550 Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
551 while(nclInnLay--){ //loop starting from layer innLay
553 Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
559 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
562 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
563 &trs,primaryVertex[2],pflag);
564 for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
566 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
567 &trs,primaryVertex[2],pflag);
571 fPoint3[0]=fPointc[0];
572 fPoint3[1]=fPointc[1];
580 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
581 if(nClusLay[nnp]!=0) layOK+=1;
584 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]));
585 AliITStrackV2* tr2 = 0;
586 Bool_t onePoint = kTRUE;
587 tr2 = FitTrack(&trs,primaryVertex,onePoint);
591 AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
593 StoreTrack(tr2,event,useAllClusters);
598 }//end loop on clusters of innLay
599 } //end loop on window sizes
601 } //end loop on innLay
602 } // end search 1-point tracks
604 if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
610 //________________________________________________________________________
612 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
613 //fit of the found track (most general case, also <6 points, layers missing)
614 // A.Dainese 16.11.07
617 const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
619 static Int_t firstmod[AliITSgeomTGeo::kNLayers];
620 static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
621 static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
622 static Int_t end[AliITSgeomTGeo::kNLayers];
623 static Int_t indices[AliITSgeomTGeo::kNLayers];
625 static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
627 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
628 firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
630 for(Int_t j=0;j<kMaxClu; j++){
638 Int_t nclusters = tr->GetNumberOfClustersSA();
639 for(Int_t ncl=0;ncl<nclusters;ncl++){
640 Int_t index = tr->GetClusterIndexSA(ncl);
641 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
642 Int_t lay = (index & 0xf0000000) >> 28;
643 Int_t nInLay=end[lay];
644 listlayer[lay][nInLay]=cl;
645 clind[lay][nInLay]=index;
649 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
650 for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
651 Int_t mark = tr->GetClusterMark(nlay,ncl);
652 clmark[nlay][ncl]=mark;
657 Int_t firstLay=-1,secondLay=-1;
658 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
664 } else if(secondLay==-1) {
670 if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
671 TClonesArray &arrMI= *fListOfTracks;
672 TClonesArray &arrSA= *fListOfSATracks;
673 Int_t nFoundTracks=0;
676 for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
678 for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
680 for(Int_t l2=0;l2<end[2];l2++){ //loop on layer 3
682 for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
684 for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
686 for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
689 // estimate curvature from 2 innermost points (or innermost point + vertex)
691 Int_t iFirstLay=indices[firstLay];
692 Int_t mrk1=clmark[firstLay][iFirstLay];
694 AliITSRecPoint* p1=(AliITSRecPoint*)listlayer[firstLay][iFirstLay];
695 Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
696 Int_t layer,ladder,detector;
697 AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
698 Float_t yclu1 = p1->GetY();
699 Float_t zclu1 = p1->GetZ();
703 Double_t cv=0,tgl2=0,phi2=0;
704 AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
710 Int_t iSecondLay=indices[secondLay];
711 Int_t mrk2=clmark[secondLay][iSecondLay];
712 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
716 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
717 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
718 phi2 = TMath::ATan2((y2-y1),(x2-x1));
719 } else { // special case of 1-point tracks, only for cosmics (B=0)
720 x2 = primaryVertex[0];
721 y2 = primaryVertex[1];
722 z2 = primaryVertex[2];
724 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
725 phi2 = TMath::ATan2((y1-y2),(x1-x2));
728 // create track and attach it the RecPoints
729 AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
730 for(Int_t iLay=5; iLay>=0; iLay--){
731 Int_t iInLay=indices[iLay];
732 AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay];
734 trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
735 trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
739 //fit with Kalman filter using AliITStrackerMI::RefitAt()
740 AliITStrackSA ot(trac);
742 ot.ResetCovariance(10.);
745 // Propagate inside the innermost layer with a cluster
746 if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
748 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6
749 AliITStrackMI otrack2(ot);
750 otrack2.ResetCovariance(10.);
751 otrack2.ResetClusters();
752 //fit from layer 6 to layer 1
753 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) {
754 new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
755 new(arrSA[nFoundTracks]) AliITStrackSA(trac);
771 if(fListOfTracks->GetEntries()==0) return 0;
773 Int_t lowchi2 = FindTrackLowChiSquare();
774 AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
775 AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
777 if(otrack==0) return 0;
779 CookLabel(otrack,0.); //MI change - to see fake ratio
780 Int_t label=FindLabel(otrack);
781 otrack->SetLabel(label);
784 otrack->CookdEdx(low,up);
786 //remove clusters of found track
787 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
788 for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
789 Int_t index = trsa->GetClusterMark(nlay,cln);
790 fCluLayer[nlay]->RemoveAt(index);
791 RemoveClusterCoord(nlay,index);
792 fCluLayer[nlay]->Compress();
800 //_______________________________________________________
801 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const
804 // Add new track to the ESD
806 AliESDtrack outtrack;
807 outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
808 if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
809 for(Int_t i=0;i<12;i++) {
810 outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
812 Double_t sdedx[4]={0.,0.,0.,0.};
813 for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
814 outtrack.SetITSdEdxSamples(sdedx);
815 event->AddTrack(&outtrack);
821 //_______________________________________________________
822 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
823 //function used to to find the clusters associated to the track
825 if(ForceSkippingOfLayer(layer)) return 0;
828 AliITSlayer &lay = fgLayers[layer];
829 Double_t r=lay.GetR();
831 Float_t cx1,cx2,cy1,cy2;
832 FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
833 if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
835 Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
836 Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
837 fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
841 Double_t phiExpect=fPhiEstimate;
842 Double_t lamExpect=fLambdac;
844 Int_t ncl = fCluLayer[layer]->GetEntriesFast();
845 for (Int_t index=0; index<ncl; index++) {
846 AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->UncheckedAt(index);
849 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
851 Double_t lambda = arr->GetLambda();
852 if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
854 Double_t phi = arr->GetPhi();
855 Double_t deltaPhi = phi-phiExpect;
856 if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
857 else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
858 if (TMath::Abs(deltaPhi)>phiwindow) continue;
860 if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
861 if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
862 Int_t orind = arr->GetOrInd();
863 trs->AddClusterSA(layer,orind);
864 trs->AddClusterMark(layer,index);
869 fPointc[0]=arr->GetX();
870 fPointc[1]=arr->GetY();
876 //________________________________________________________________
877 Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
878 // Sets the first point (seed) for tracking
880 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[lay]->UncheckedAt(clu);
881 if(!cl) return kFALSE;
882 if (cl->GetQ()<=0) return kFALSE;
883 if(lay>1 && cl->GetQ()<=fMinQ) return kFALSE;
885 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
886 fPhic = arr->GetPhi();
887 fLambdac = arr->GetLambda();
888 fPhiEstimate = fPhic;
889 fPoint1[0]=primaryVertex[0];
890 fPoint1[1]=primaryVertex[1];
891 fPoint2[0]=arr->GetX();
892 fPoint2[1]=arr->GetY();
896 //________________________________________________________________
897 void AliITStrackerSA::UpdatePoints(){
898 //update of points for the estimation of the curvature
900 fPoint2[0]=fPoint3[0];
901 fPoint2[1]=fPoint3[1];
902 fPoint3[0]=fPointc[0];
903 fPoint3[1]=fPointc[1];
908 //___________________________________________________________________
909 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){
911 //given (x,y) of three recpoints (in global coordinates)
912 //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
914 Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
916 a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
917 b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
918 c = -x1*x1-y1*y1-a*x1-b*y1;
921 //__________________________________________________________________________
922 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){
924 //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
925 //c2 is -rlayer*rlayer
929 Double_t aA = (b1*b1)/(a1*a1)+1;
930 Double_t bB = (-2*m*b1/(a1*a1));
931 Double_t cC = c2+(m*m)/(a1*a1);
932 Double_t dD = bB*bB-4*aA*cC;
935 y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
936 y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
937 x1 = (c2-c1-b1*y1)/a1;
938 x2 = (c2-c1-b1*y2)/a1;
942 //____________________________________________________________________
943 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
944 x2,Double_t y2,Double_t x3,Double_t y3){
946 //calculates the curvature of track
947 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
949 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
950 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
951 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
954 if((a*a+b*b-4*c)<0) return 0;
955 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
958 if((x1>0 && y1>0 && x1<xc)) rad*=-1;
959 if((x1<0 && y1>0 && x1<xc)) rad*=-1;
960 // if((x1<0 && y1<0 && x1<xc)) rad*=-1;
961 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
968 //____________________________________________________________________
969 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
971 //Returns the point closest to pp
973 Double_t diff1 = p1-pp;
974 Double_t diff2 = p2-pp;
976 if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
977 else fPhiEstimate=p2;
983 //_________________________________________________________________
984 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
985 // returns track with lowest chi square
986 Int_t dim=fListOfTracks->GetEntries();
988 AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
989 Double_t minChi2=trk->GetChi2();
991 for(Int_t i=1;i<dim;i++){
992 trk = (AliITStrackV2*)fListOfTracks->At(i);
993 Double_t chi2=trk->GetChi2();
1002 //__________________________________________________________
1003 Int_t AliITStrackerSA::FindLabel(AliITStrackV2* track){
1006 Int_t labl[AliITSgeomTGeo::kNLayers][3];
1007 Int_t cnts[AliITSgeomTGeo::kNLayers][3];
1008 for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
1009 for(Int_t k=0;k<3;k++){
1015 for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
1016 Int_t indexc = track->GetClusterIndex(i);
1017 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
1018 Int_t iLayer=cl->GetLayer();
1019 for(Int_t k=0;k<3;k++){
1020 labl[iLayer][k]=cl->GetLabel(k);
1021 if(labl[iLayer][k]<0) iNotLabel++;
1024 if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
1026 for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
1027 for(Int_t j2=0; j2<j1; j2++){
1028 for(Int_t k1=0; k1<3; k1++){
1029 for(Int_t k2=0; k2<3; k2++){
1030 if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
1042 for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
1043 for(Int_t k=0;k<3;k++){
1044 if(cnts[j][k]>cntMax && labl[j][k]>=0){
1052 for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
1053 if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
1055 if(lflag<track->GetNumberOfClusters()) label = -label;
1058 //_____________________________________________________________________________
1059 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
1060 // Set sizes of the phi and lambda windows used for track finding
1062 if(fPhiWin) delete [] fPhiWin;
1063 if(fLambdaWin) delete [] fLambdaWin;
1064 fPhiWin = new Double_t[fNloop];
1065 fLambdaWin = new Double_t[fNloop];
1066 Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
1067 Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
1068 for(Int_t k=0;k<fNloop;k++){
1069 Float_t phi=phimin+k*stepPhi;
1070 Float_t lam=lambdamin+k*stepLambda;
1075 //_____________________________________________________________________________
1076 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1077 // Set sizes of the phi and lambda windows used for track finding
1079 if(phi){ // user defined values
1080 fPhiWin = new Double_t[fNloop];
1081 fLambdaWin = new Double_t[fNloop];
1082 for(Int_t k=0;k<fNloop;k++){
1084 fLambdaWin[k]=lam[k];
1087 else { // default values
1089 Double_t phid[33] = {0.002,0.003,0.004,0.0045,0.0047,
1090 0.005,0.0053,0.0055,
1091 0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1092 0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1093 0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1094 Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1096 0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1097 0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1098 0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1105 fPhiWin = new Double_t[fNloop];
1106 fLambdaWin = new Double_t[fNloop];
1108 Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1110 for(Int_t k=0;k<fNloop;k++){
1111 fPhiWin[k]=phid[k]*factor;
1112 fLambdaWin[k]=lambdad[k]*factor;
1118 //_______________________________________________________________________
1119 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1120 //Returns values of phi (azimuthal) and lambda angles for a given cluster
1122 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1123 Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1124 Float_t tx,ty,tz; fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1126 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1127 Double_t phi1=TMath::Pi()/2+alpha;
1128 if (lay==1) phi1+=TMath::Pi();
1130 Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1131 Float_t r=tx*cp+ty*sp;
1133 xyz= r*cp - cl->GetY()*sp;
1134 y= r*sp + cl->GetY()*cp;
1138 cl->GetGlobalXYZ(xyz);
1143 phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1144 lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1147 //________________________________________________________________________
1148 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1150 //returns sigmax, y, z of cluster in global coordinates
1152 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1154 AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1156 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1157 Double_t phi=TMath::Pi()/2+alpha;
1158 if (lay==1) phi+=TMath::Pi();
1160 Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1163 cl->GetGlobalCov(covm);
1164 sx=TMath::Sqrt(covm[0]);
1165 sy=TMath::Sqrt(covm[3]);
1166 sz=TMath::Sqrt(covm[5]);
1168 sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1169 sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1170 sz = TMath::Sqrt(cl->GetSigmaZ2());