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;
494 if(nextLay<0 || nextLay == 6) goon = kFALSE;
497 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
498 &trs,primaryVertex[2],pflag);
502 fPoint3[0]=fPointc[0];
503 fPoint3[1]=fPointc[1];
509 if(nextLay<0 || nextLay==6) goon=kFALSE;
515 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-theLay;nnp++){
516 if(nClusLay[nnp]!=0) layOK+=1;
519 for(Int_t nnp=theLay; nnp>=0; nnp--){
520 if(nClusLay[nnp]!=0) layOK+=1;
523 if(layOK>=minNPoints){
524 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]));
525 AliITStrackV2* tr2 = 0;
526 tr2 = FitTrack(&trs,primaryVertex);
530 AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
532 StoreTrack(tr2,event,useAllClusters);
537 }//end loop on clusters of theLay
538 } //end loop on window sizes
539 } //end loop on theLay
540 }//end loop on min points
542 // search for 1-point tracks in SPD, only for cosmics
543 // (A.Dainese 21.03.08)
544 if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
545 TMath::Abs(event->GetMagneticField())<0.01) {
546 Int_t outerLayer=1; // only SPD
547 for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
548 // counter for clusters on each layer
550 for(Int_t nloop=0;nloop<fNloop;nloop++){
551 Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
552 while(nclInnLay--){ //loop starting from layer innLay
554 Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
560 for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
563 nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
564 &trs,primaryVertex[2],pflag);
565 for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
567 nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
568 &trs,primaryVertex[2],pflag);
572 fPoint3[0]=fPointc[0];
573 fPoint3[1]=fPointc[1];
581 for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
582 if(nClusLay[nnp]!=0) layOK+=1;
585 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]));
586 AliITStrackV2* tr2 = 0;
587 Bool_t onePoint = kTRUE;
588 tr2 = FitTrack(&trs,primaryVertex,onePoint);
592 AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
594 StoreTrack(tr2,event,useAllClusters);
599 }//end loop on clusters of innLay
600 } //end loop on window sizes
602 } //end loop on innLay
603 } // end search 1-point tracks
605 if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
611 //________________________________________________________________________
613 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
614 //fit of the found track (most general case, also <6 points, layers missing)
615 // A.Dainese 16.11.07
618 const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
620 static Int_t firstmod[AliITSgeomTGeo::kNLayers];
621 static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
622 static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
623 static Int_t end[AliITSgeomTGeo::kNLayers];
624 static Int_t indices[AliITSgeomTGeo::kNLayers];
626 static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
628 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
629 firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
631 for(Int_t j=0;j<kMaxClu; j++){
639 Int_t nclusters = tr->GetNumberOfClustersSA();
640 for(Int_t ncl=0;ncl<nclusters;ncl++){
641 Int_t index = tr->GetClusterIndexSA(ncl);
642 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
643 Int_t lay = (index & 0xf0000000) >> 28;
644 Int_t nInLay=end[lay];
645 listlayer[lay][nInLay]=cl;
646 clind[lay][nInLay]=index;
650 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
651 for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
652 Int_t mark = tr->GetClusterMark(nlay,ncl);
653 clmark[nlay][ncl]=mark;
658 Int_t firstLay=-1,secondLay=-1;
659 for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
665 } else if(secondLay==-1) {
671 if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
672 TClonesArray &arrMI= *fListOfTracks;
673 TClonesArray &arrSA= *fListOfSATracks;
674 Int_t nFoundTracks=0;
677 for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
679 for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
681 for(Int_t l2=0;l2<end[2];l2++){ //loop on layer 3
683 for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
685 for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
687 for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
690 // estimate curvature from 2 innermost points (or innermost point + vertex)
692 Int_t iFirstLay=indices[firstLay];
693 Int_t mrk1=clmark[firstLay][iFirstLay];
695 AliITSRecPoint* p1=(AliITSRecPoint*)listlayer[firstLay][iFirstLay];
696 Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
697 Int_t layer,ladder,detector;
698 AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
699 Float_t yclu1 = p1->GetY();
700 Float_t zclu1 = p1->GetZ();
704 Double_t cv=0,tgl2=0,phi2=0;
705 AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
711 Int_t iSecondLay=indices[secondLay];
712 Int_t mrk2=clmark[secondLay][iSecondLay];
713 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
717 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
718 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
719 phi2 = TMath::ATan2((y2-y1),(x2-x1));
720 } else { // special case of 1-point tracks, only for cosmics (B=0)
721 x2 = primaryVertex[0];
722 y2 = primaryVertex[1];
723 z2 = primaryVertex[2];
725 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
726 phi2 = TMath::ATan2((y1-y2),(x1-x2));
729 // create track and attach it the RecPoints
730 AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
731 for(Int_t iLay=5; iLay>=0; iLay--){
732 Int_t iInLay=indices[iLay];
733 AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay];
735 trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
736 trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
740 //fit with Kalman filter using AliITStrackerMI::RefitAt()
741 AliITStrackSA ot(trac);
743 ot.ResetCovariance(10.);
746 // Propagate inside the innermost layer with a cluster
747 if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
749 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6
750 AliITStrackMI otrack2(ot);
751 otrack2.ResetCovariance(10.);
752 otrack2.ResetClusters();
753 //fit from layer 6 to layer 1
754 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) {
755 new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
756 new(arrSA[nFoundTracks]) AliITStrackSA(trac);
772 if(fListOfTracks->GetEntries()==0) return 0;
774 Int_t lowchi2 = FindTrackLowChiSquare();
775 AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
776 AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
778 if(otrack==0) return 0;
780 CookLabel(otrack,0.); //MI change - to see fake ratio
781 Int_t label=FindLabel(otrack);
782 otrack->SetLabel(label);
785 otrack->CookdEdx(low,up);
787 //remove clusters of found track
788 for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
789 for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
790 Int_t index = trsa->GetClusterMark(nlay,cln);
791 fCluLayer[nlay]->RemoveAt(index);
792 RemoveClusterCoord(nlay,index);
793 fCluLayer[nlay]->Compress();
801 //_______________________________________________________
802 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const
805 // Add new track to the ESD
807 AliESDtrack outtrack;
808 outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
809 if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
810 for(Int_t i=0;i<12;i++) {
811 outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
813 Double_t sdedx[4]={0.,0.,0.,0.};
814 for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
815 outtrack.SetITSdEdxSamples(sdedx);
816 event->AddTrack(&outtrack);
822 //_______________________________________________________
823 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
824 //function used to to find the clusters associated to the track
826 if(ForceSkippingOfLayer(layer)) return 0;
829 AliITSlayer &lay = fgLayers[layer];
830 Double_t r=lay.GetR();
832 Float_t cx1,cx2,cy1,cy2;
833 FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
834 if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
836 Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
837 Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
838 fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
842 Double_t phiExpect=fPhiEstimate;
843 Double_t lamExpect=fLambdac;
845 Int_t ncl = fCluLayer[layer]->GetEntriesFast();
846 for (Int_t index=0; index<ncl; index++) {
847 AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->UncheckedAt(index);
850 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
852 Double_t lambda = arr->GetLambda();
853 if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
855 Double_t phi = arr->GetPhi();
856 Double_t deltaPhi = phi-phiExpect;
857 if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
858 else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
859 if (TMath::Abs(deltaPhi)>phiwindow) continue;
861 if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
862 if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
863 Int_t orind = arr->GetOrInd();
864 trs->AddClusterSA(layer,orind);
865 trs->AddClusterMark(layer,index);
870 fPointc[0]=arr->GetX();
871 fPointc[1]=arr->GetY();
877 //________________________________________________________________
878 Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
879 // Sets the first point (seed) for tracking
881 AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[lay]->UncheckedAt(clu);
882 if(!cl) return kFALSE;
883 if (cl->GetQ()<=0) return kFALSE;
884 if(lay>1 && cl->GetQ()<=fMinQ) return kFALSE;
886 AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
887 fPhic = arr->GetPhi();
888 fLambdac = arr->GetLambda();
889 fPhiEstimate = fPhic;
890 fPoint1[0]=primaryVertex[0];
891 fPoint1[1]=primaryVertex[1];
892 fPoint2[0]=arr->GetX();
893 fPoint2[1]=arr->GetY();
897 //________________________________________________________________
898 void AliITStrackerSA::UpdatePoints(){
899 //update of points for the estimation of the curvature
901 fPoint2[0]=fPoint3[0];
902 fPoint2[1]=fPoint3[1];
903 fPoint3[0]=fPointc[0];
904 fPoint3[1]=fPointc[1];
909 //___________________________________________________________________
910 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){
912 //given (x,y) of three recpoints (in global coordinates)
913 //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
915 Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
917 a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
918 b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
919 c = -x1*x1-y1*y1-a*x1-b*y1;
922 //__________________________________________________________________________
923 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){
925 //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
926 //c2 is -rlayer*rlayer
930 Double_t aA = (b1*b1)/(a1*a1)+1;
931 Double_t bB = (-2*m*b1/(a1*a1));
932 Double_t cC = c2+(m*m)/(a1*a1);
933 Double_t dD = bB*bB-4*aA*cC;
936 y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
937 y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
938 x1 = (c2-c1-b1*y1)/a1;
939 x2 = (c2-c1-b1*y2)/a1;
943 //____________________________________________________________________
944 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
945 x2,Double_t y2,Double_t x3,Double_t y3){
947 //calculates the curvature of track
948 Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
950 Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
951 Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
952 Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
955 if((a*a+b*b-4*c)<0) return 0;
956 Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
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;
962 // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
969 //____________________________________________________________________
970 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
972 //Returns the point closest to pp
974 Double_t diff1 = p1-pp;
975 Double_t diff2 = p2-pp;
977 if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
978 else fPhiEstimate=p2;
984 //_________________________________________________________________
985 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
986 // returns track with lowest chi square
987 Int_t dim=fListOfTracks->GetEntries();
989 AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
990 Double_t minChi2=trk->GetChi2();
992 for(Int_t i=1;i<dim;i++){
993 trk = (AliITStrackV2*)fListOfTracks->At(i);
994 Double_t chi2=trk->GetChi2();
1003 //__________________________________________________________
1004 Int_t AliITStrackerSA::FindLabel(AliITStrackV2* track){
1007 Int_t labl[AliITSgeomTGeo::kNLayers][3];
1008 Int_t cnts[AliITSgeomTGeo::kNLayers][3];
1009 for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
1010 for(Int_t k=0;k<3;k++){
1016 for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
1017 Int_t indexc = track->GetClusterIndex(i);
1018 AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
1019 Int_t iLayer=cl->GetLayer();
1020 for(Int_t k=0;k<3;k++){
1021 labl[iLayer][k]=cl->GetLabel(k);
1022 if(labl[iLayer][k]<0) iNotLabel++;
1025 if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
1027 for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
1028 for(Int_t j2=0; j2<j1; j2++){
1029 for(Int_t k1=0; k1<3; k1++){
1030 for(Int_t k2=0; k2<3; k2++){
1031 if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
1043 for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
1044 for(Int_t k=0;k<3;k++){
1045 if(cnts[j][k]>cntMax && labl[j][k]>=0){
1053 for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
1054 if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
1056 if(lflag<track->GetNumberOfClusters()) label = -label;
1059 //_____________________________________________________________________________
1060 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
1061 // Set sizes of the phi and lambda windows used for track finding
1063 if(fPhiWin) delete [] fPhiWin;
1064 if(fLambdaWin) delete [] fLambdaWin;
1065 fPhiWin = new Double_t[fNloop];
1066 fLambdaWin = new Double_t[fNloop];
1067 Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
1068 Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
1069 for(Int_t k=0;k<fNloop;k++){
1070 Float_t phi=phimin+k*stepPhi;
1071 Float_t lam=lambdamin+k*stepLambda;
1076 //_____________________________________________________________________________
1077 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1078 // Set sizes of the phi and lambda windows used for track finding
1080 if(phi){ // user defined values
1081 fPhiWin = new Double_t[fNloop];
1082 fLambdaWin = new Double_t[fNloop];
1083 for(Int_t k=0;k<fNloop;k++){
1085 fLambdaWin[k]=lam[k];
1088 else { // default values
1090 Double_t phid[33] = {0.002,0.003,0.004,0.0045,0.0047,
1091 0.005,0.0053,0.0055,
1092 0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1093 0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1094 0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1095 Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1097 0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1098 0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1099 0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1106 fPhiWin = new Double_t[fNloop];
1107 fLambdaWin = new Double_t[fNloop];
1109 Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1111 for(Int_t k=0;k<fNloop;k++){
1112 fPhiWin[k]=phid[k]*factor;
1113 fLambdaWin[k]=lambdad[k]*factor;
1119 //_______________________________________________________________________
1120 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1121 //Returns values of phi (azimuthal) and lambda angles for a given cluster
1123 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1124 Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1125 Float_t tx,ty,tz; fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1127 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1128 Double_t phi1=TMath::Pi()/2+alpha;
1129 if (lay==1) phi1+=TMath::Pi();
1131 Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1132 Float_t r=tx*cp+ty*sp;
1134 xyz= r*cp - cl->GetY()*sp;
1135 y= r*sp + cl->GetY()*cp;
1139 cl->GetGlobalXYZ(xyz);
1144 phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1145 lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1148 //________________________________________________________________________
1149 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1151 //returns sigmax, y, z of cluster in global coordinates
1153 Double_t rot[9]; fGeom->GetRotMatrix(module,rot);
1155 AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1157 Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1158 Double_t phi=TMath::Pi()/2+alpha;
1159 if (lay==1) phi+=TMath::Pi();
1161 Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1164 cl->GetGlobalCov(covm);
1165 sx=TMath::Sqrt(covm[0]);
1166 sy=TMath::Sqrt(covm[3]);
1167 sz=TMath::Sqrt(covm[5]);
1169 sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1170 sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1171 sz = TMath::Sqrt(cl->GetSigmaZ2());