]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITStrackerMI.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity warnings
[u/mrichter/AliRoot.git] / ITS / AliITStrackerMI.cxx
CommitLineData
e43c066c 1/**************************************************************************
9a5326b1 2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
e43c066c 3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
fddf8459 15/* $Id$ */
bf6adc12 16
e43c066c 17//-------------------------------------------------------------------------
18// Implementation of the ITS tracker class
00a7cc50 19// It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
e43c066c 20// and fills with them the ESD
15dd636f 21// Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
23197852 22// Current support and development:
23// Andrea Dainese, andrea.dainese@lnl.infn.it
afd25725 24// dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25// Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26// Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
e43c066c 27//-------------------------------------------------------------------------
bf6adc12 28
29#include <TMatrixD.h>
e43c066c 30#include <TTree.h>
afd25725 31#include <TDatabasePDG.h>
9a5326b1 32#include <TString.h>
e50912db 33#include <TRandom.h>
a70ed6ad 34#include <TTreeStream.h>
879cdb02 35#include <TVector3.h>
bf6adc12 36
4187a601 37#include "AliLog.h"
88e3c173 38#include "AliGeomManager.h"
a70ed6ad 39#include "AliITSPlaneEff.h"
a70ed6ad 40#include "AliTrackPointArray.h"
af885e0f 41#include "AliESDEvent.h"
aad72f45 42#include "AliESDtrack.h"
6c94f330 43#include "AliV0.h"
a70ed6ad 44#include "AliITSChannelStatus.h"
45#include "AliITSDetTypeRec.h"
00a7cc50 46#include "AliITSRecPoint.h"
2b295109 47#include "AliITSRecPointContainer.h"
e341247d 48#include "AliITSgeomTGeo.h"
44347160 49#include "AliITSReconstructor.h"
572f41f9 50#include "AliITSClusterParam.h"
23197852 51#include "AliITSsegmentation.h"
52#include "AliITSCalibration.h"
7167ae53 53#include "AliITSPlaneEffSPD.h"
54#include "AliITSPlaneEffSDD.h"
55#include "AliITSPlaneEffSSD.h"
44396017 56#include "AliITSV0Finder.h"
ae00569a 57#include "AliITStrackerMI.h"
bb7e41dd 58#include "AliMathBase.h"
e43c066c 59
60ClassImp(AliITStrackerMI)
e43c066c 61
e50912db 62AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
44347160 63
8221b41b 64AliITStrackerMI::AliITStrackerMI():AliTracker(),
65fI(0),
66fBestTrack(),
67fTrackToFollow(),
68fTrackHypothesys(),
69fBestHypothesys(),
70fOriginal(),
71fCurrentEsdTrack(),
72fPass(0),
73fAfterV0(kFALSE),
74fLastLayerToTrackTo(0),
afd25725 75fCoefficients(0),
8221b41b 76fEsd(0),
e50912db 77fTrackingPhase("Default"),
6518a6c5 78fUseTGeo(3),
e50912db 79fNtracks(0),
80fxOverX0Pipe(-1.),
81fxTimesRhoPipe(-1.),
82fxOverX0PipeTrks(0),
83fxTimesRhoPipeTrks(0),
84fxOverX0ShieldTrks(0),
85fxTimesRhoShieldTrks(0),
86fxOverX0LayerTrks(0),
87fxTimesRhoLayerTrks(0),
4a66240a 88fDebugStreamer(0),
23197852 89fITSChannelStatus(0),
a70ed6ad 90fkDetTypeRec(0),
ae00569a 91fPlaneEff(0) {
8221b41b 92 //Default constructor
e50912db 93 Int_t i;
94 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
225cff59 95 for(i=0;i<2;i++) {
96 fxOverX0Shield[i]=-1.;
97 fxTimesRhoShield[i]=-1.;
98 fConstraint[i]=0;
99 }
e50912db 100 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
a0b53cc4 101 fOriginal.SetOwner();
225cff59 102 for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
103 for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
104
8221b41b 105}
44347160 106//------------------------------------------------------------------------
1f3e997f 107AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
e50912db 108fI(AliITSgeomTGeo::GetNLayers()),
8221b41b 109fBestTrack(),
110fTrackToFollow(),
111fTrackHypothesys(),
112fBestHypothesys(),
113fOriginal(),
114fCurrentEsdTrack(),
115fPass(0),
116fAfterV0(kFALSE),
e50912db 117fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
afd25725 118fCoefficients(0),
8221b41b 119fEsd(0),
e50912db 120fTrackingPhase("Default"),
6518a6c5 121fUseTGeo(3),
e50912db 122fNtracks(0),
123fxOverX0Pipe(-1.),
124fxTimesRhoPipe(-1.),
125fxOverX0PipeTrks(0),
126fxTimesRhoPipeTrks(0),
127fxOverX0ShieldTrks(0),
128fxTimesRhoShieldTrks(0),
129fxOverX0LayerTrks(0),
130fxTimesRhoLayerTrks(0),
4a66240a 131fDebugStreamer(0),
23197852 132fITSChannelStatus(0),
a70ed6ad 133fkDetTypeRec(0),
ae00569a 134fPlaneEff(0) {
e43c066c 135 //--------------------------------------------------------------------
136 //This is the AliITStrackerMI constructor
137 //--------------------------------------------------------------------
e341247d 138 if (geom) {
139 AliWarning("\"geom\" is actually a dummy argument !");
140 }
141
a0b53cc4 142 fOriginal.SetOwner();
afd25725 143 fCoefficients = 0;
628e7bb0 144 fAfterV0 = kFALSE;
e341247d 145
e50912db 146 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
e341247d 147 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
148 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
149
150 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
bf1cd124 151 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
e43c066c 152 Double_t poff=TMath::ATan2(y,x);
153 Double_t zoff=z;
bf1cd124 154
155 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
e341247d 156 Double_t r=TMath::Sqrt(x*x + y*y);
e43c066c 157
e341247d 158 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
e43c066c 159 r += TMath::Sqrt(x*x + y*y);
e341247d 160 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
e43c066c 161 r += TMath::Sqrt(x*x + y*y);
e341247d 162 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
e43c066c 163 r += TMath::Sqrt(x*x + y*y);
164 r*=0.25;
165
166 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
167
168 for (Int_t j=1; j<nlad+1; j++) {
169 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
e341247d 170 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
1f3e997f 171 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
172 m.Multiply(tm);
a5d0fae6 173 Double_t txyz[3]={0.};
174 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
1f3e997f 175 m.LocalToMaster(txyz,xyz);
a5d0fae6 176 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
1f3e997f 177 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
e341247d 178
179 if (phi<0) phi+=TMath::TwoPi();
180 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
181
e43c066c 182 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
183 new(&det) AliITSdetector(r,phi);
1c97ce2f 184 // compute the real radius (with misalignment)
185 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
186 mmisal.Multiply(tm);
187 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
188 mmisal.LocalToMaster(txyz,xyz);
189 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
190 det.SetRmisal(rmisal);
191
23197852 192 } // end loop on detectors
193 } // end loop on ladders
9cd19c63 194 fForceSkippingOfLayer[i-1] = 0;
23197852 195 } // end loop on layers
e43c066c 196
e43c066c 197
e50912db 198 fI=AliITSgeomTGeo::GetNLayers();
e43c066c 199
200 fPass=0;
201 fConstraint[0]=1; fConstraint[1]=0;
202
afd25725 203 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
204 AliITSReconstructor::GetRecoParam()->GetYVdef(),
205 AliITSReconstructor::GetRecoParam()->GetZVdef()};
206 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
207 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
208 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
209 SetVertex(xyzVtx,ersVtx);
e43c066c 210
e50912db 211 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
628e7bb0 212 for (Int_t i=0;i<100000;i++){
213 fBestTrackIndex[i]=0;
214 }
afd25725 215
216 // store positions of centre of SPD modules (in z)
bf1cd124 217 // The convetion is that fSPDdetzcentre is ordered from -z to +z
afd25725 218 Double_t tr[3];
c3fd826e 219 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
bf1cd124 220 if (tr[2]<0) { // old geom (up to v5asymmPPR)
221 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
222 fSPDdetzcentre[0] = tr[2];
223 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
224 fSPDdetzcentre[1] = tr[2];
225 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
226 fSPDdetzcentre[2] = tr[2];
227 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
228 fSPDdetzcentre[3] = tr[2];
229 } else { // new geom (from v11Hybrid)
230 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
231 fSPDdetzcentre[0] = tr[2];
232 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
233 fSPDdetzcentre[1] = tr[2];
234 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
235 fSPDdetzcentre[2] = tr[2];
236 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
237 fSPDdetzcentre[3] = tr[2];
238 }
afd25725 239
afd25725 240 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
6518a6c5 241 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
242 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
243 fUseTGeo = 3;
244 }
e50912db 245
246 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
247 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
248
c0b60eb0 249 if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
250 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
81e97e0d 251
7167ae53 252 // only for plane efficiency evaluation
58e8dc31 253 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
254 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
0ed58a47 255 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
1a64bda6 256 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
0ed58a47 257 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
258 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
259 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
260 else fPlaneEff = new AliITSPlaneEffSSD();
275a301c 261 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
0ed58a47 262 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
263 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
7167ae53 264 }
628e7bb0 265}
225cff59 266/*
44347160 267//------------------------------------------------------------------------
8221b41b 268AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
269fI(tracker.fI),
270fBestTrack(tracker.fBestTrack),
271fTrackToFollow(tracker.fTrackToFollow),
272fTrackHypothesys(tracker.fTrackHypothesys),
273fBestHypothesys(tracker.fBestHypothesys),
274fOriginal(tracker.fOriginal),
275fCurrentEsdTrack(tracker.fCurrentEsdTrack),
276fPass(tracker.fPass),
277fAfterV0(tracker.fAfterV0),
278fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
afd25725 279fCoefficients(tracker.fCoefficients),
8221b41b 280fEsd(tracker.fEsd),
e50912db 281fTrackingPhase(tracker.fTrackingPhase),
afd25725 282fUseTGeo(tracker.fUseTGeo),
e50912db 283fNtracks(tracker.fNtracks),
284fxOverX0Pipe(tracker.fxOverX0Pipe),
285fxTimesRhoPipe(tracker.fxTimesRhoPipe),
286fxOverX0PipeTrks(0),
287fxTimesRhoPipeTrks(0),
288fxOverX0ShieldTrks(0),
289fxTimesRhoShieldTrks(0),
290fxOverX0LayerTrks(0),
291fxTimesRhoLayerTrks(0),
4a66240a 292fDebugStreamer(tracker.fDebugStreamer),
23197852 293fITSChannelStatus(tracker.fITSChannelStatus),
a70ed6ad 294fkDetTypeRec(tracker.fkDetTypeRec),
ae00569a 295fPlaneEff(tracker.fPlaneEff) {
8221b41b 296 //Copy constructor
a0b53cc4 297 fOriginal.SetOwner();
e50912db 298 Int_t i;
299 for(i=0;i<4;i++) {
300 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
301 }
302 for(i=0;i<6;i++) {
303 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
304 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
305 }
306 for(i=0;i<2;i++) {
307 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
308 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
309 }
8221b41b 310}
225cff59 311
44347160 312//------------------------------------------------------------------------
8221b41b 313AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
314 //Assignment operator
315 this->~AliITStrackerMI();
316 new(this) AliITStrackerMI(tracker);
317 return *this;
318}
225cff59 319*/
44347160 320//------------------------------------------------------------------------
628e7bb0 321AliITStrackerMI::~AliITStrackerMI()
322{
323 //
324 //destructor
325 //
e50912db 326 if (fCoefficients) delete [] fCoefficients;
327 DeleteTrksMaterialLUT();
81e97e0d 328 if (fDebugStreamer) {
329 //fDebugStreamer->Close();
330 delete fDebugStreamer;
331 }
23197852 332 if(fITSChannelStatus) delete fITSChannelStatus;
333 if(fPlaneEff) delete fPlaneEff;
e43c066c 334}
44347160 335//------------------------------------------------------------------------
23197852 336void AliITStrackerMI::ReadBadFromDetTypeRec() {
337 //--------------------------------------------------------------------
338 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
339 //i.e. from OCDB
340 //--------------------------------------------------------------------
341
342 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
343
4187a601 344 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
23197852 345
a70ed6ad 346 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
23197852 347
348 // ITS channels map
349 if(fITSChannelStatus) delete fITSChannelStatus;
a70ed6ad 350 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
23197852 351
352 // ITS detectors and chips
353 Int_t i=0,j=0,k=0,ndet=0;
354 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
4187a601 355 Int_t nBadDetsPerLayer=0;
23197852 356 ndet=AliITSgeomTGeo::GetNDetectors(i);
357 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
358 for (k=1; k<ndet+1; k++) {
359 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
a70ed6ad 360 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
4187a601 361 if(det.IsBad()) {nBadDetsPerLayer++;}
23197852 362 } // end loop on detectors
363 } // end loop on ladders
e9b15b0c 364 AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
23197852 365 } // end loop on layers
366
367 return;
368}
369//------------------------------------------------------------------------
e43c066c 370Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
371 //--------------------------------------------------------------------
372 //This function loads ITS clusters
373 //--------------------------------------------------------------------
2b295109 374
375 TClonesArray *clusters = NULL;
376 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
377 clusters=rpcont->FetchClusters(0,cTree);
9cd19c63 378 if(!clusters) return 1;
379
2b295109 380 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
381 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
382 return 1;
e43c066c 383 }
23197852 384 Int_t i=0,j=0,ndet=0;
e43c066c 385 Int_t detector=0;
23197852 386 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
387 ndet=fgLayers[i].GetNdetectors();
e43c066c 388 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
389 for (; j<jmax; j++) {
2b295109 390 // if (!cTree->GetEvent(j)) continue;
391 clusters = rpcont->UncheckedGetClusters(j);
392 if(!clusters)continue;
e43c066c 393 Int_t ncl=clusters->GetEntriesFast();
394 SignDeltas(clusters,GetZ());
1f3e997f 395
e43c066c 396 while (ncl--) {
00a7cc50 397 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
75fb37cc 398 detector=c->GetDetectorIndex();
a504d56f 399
75fb37cc 400 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
6262dd3d 401
402 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
403 if(retval) {
404 AliWarning(Form("Too many clusters on layer %d!",i));
405 break;
406 }
e43c066c 407 }
6262dd3d 408
afd25725 409 // add dead zone "virtual" cluster in SPD, if there is a cluster within
410 // zwindow cm from the dead zone
bf1cd124 411 // This method assumes that fSPDdetzcentre is ordered from -z to +z
afd25725 412 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
e50912db 413 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
afd25725 414 Int_t lab[4] = {0,0,0,detector};
415 Int_t info[3] = {0,0,i};
6518a6c5 416 Float_t q = 0.; // this identifies virtual clusters
9cd19c63 417 Float_t hit[6] = {xdead,
afd25725 418 0.,
419 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
420 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
9cd19c63 421 q,
422 0.};
afd25725 423 Bool_t local = kTRUE;
424 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
e50912db 425 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 426 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
427 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
e50912db 428 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 429 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
430 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
e50912db 431 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 432 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
433 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
e50912db 434 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 435 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
436 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
e50912db 437 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 438 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
439 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
e50912db 440 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 441 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
442 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
e43c066c 443 }
6518a6c5 444 } // "virtual" clusters in SPD
e43c066c 445
446 }
447 //
448 fgLayers[i].ResetRoad(); //road defined by the cluster density
449 fgLayers[i].SortClusters();
450 }
451
4fa7f7d1 452 // check whether we have to skip some layers
453 SetForceSkippingOfLayer();
454
e43c066c 455 return 0;
456}
44347160 457//------------------------------------------------------------------------
e43c066c 458void AliITStrackerMI::UnloadClusters() {
459 //--------------------------------------------------------------------
460 //This function unloads ITS clusters
461 //--------------------------------------------------------------------
e50912db 462 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
e43c066c 463}
44347160 464//------------------------------------------------------------------------
23197852 465void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
466 //--------------------------------------------------------------------
467 // Publishes all pointers to clusters known to the tracker into the
468 // passed object array.
469 // The ownership is not transfered - the caller is not expected to delete
470 // the clusters.
471 //--------------------------------------------------------------------
472
473 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
474 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
475 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
476 array->AddLast(cl);
477 }
478 }
479
480 return;
481}
482//------------------------------------------------------------------------
44396017 483Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
e43c066c 484 //--------------------------------------------------------------------
485 // Correction for the material between the TPC and the ITS
e43c066c 486 //--------------------------------------------------------------------
e50912db 487 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
488 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
489 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
490 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
491 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
492 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
493 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
494 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
e43c066c 495 } else {
44396017 496 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
e50912db 497 return 0;
e43c066c 498 }
499
e50912db 500 return 1;
e43c066c 501}
44347160 502//------------------------------------------------------------------------
af885e0f 503Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
e43c066c 504 //--------------------------------------------------------------------
505 // This functions reconstructs ITS tracks
506 // The clusters must be already loaded !
507 //--------------------------------------------------------------------
fddf8459 508
4fa7f7d1 509 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
fddf8459 510
e50912db 511 fTrackingPhase="Clusters2Tracks";
512
e43c066c 513 TObjArray itsTracks(15000);
628e7bb0 514 fOriginal.Clear();
81e97e0d 515 fEsd = event; // store pointer to the esd
afd25725 516
9a5326b1 517 // temporary (for cosmics)
518 if(event->GetVertex()) {
519 TString title = event->GetVertex()->GetTitle();
520 if(title.Contains("cosmics")) {
521 Double_t xyz[3]={GetX(),GetY(),GetZ()};
522 Double_t exyz[3]={0.1,0.1,0.1};
523 SetVertex(xyz,exyz);
524 }
525 }
526 // temporary
75a74b11 527 Int_t noesd = 0;
e43c066c 528 {/* Read ESD tracks */
afd25725 529 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
e43c066c 530 Int_t nentr=event->GetNumberOfTracks();
75a74b11 531 noesd=nentr;
532 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
e43c066c 533 while (nentr--) {
534 AliESDtrack *esd=event->GetTrack(nentr);
4187a601 535 // ---- for debugging:
536 //if(TMath::Abs(esd->GetX()-83.65)<0.1) { FILE *f=fopen("tpc.dat","a"); fprintf(f,"%f %f %f %f %f %f\n",(Float_t)event->GetEventNumberInFile(),(Float_t)TMath::Abs(esd->GetLabel()),(Float_t)esd->GetX(),(Float_t)esd->GetY(),(Float_t)esd->GetZ(),(Float_t)esd->Pt()); fclose(f); }
e43c066c 537
538 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
539 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
540 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
51ad6848 541 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
763bdf3e 542 AliITStrackMI *t = new AliITStrackMI(*esd);
b9671574 543 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
544 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
6518a6c5 545
546
afd25725 547 // look at the ESD mass hypothesys !
548 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
e43c066c 549 // write expected q
b9671574 550 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
e43c066c 551
afd25725 552 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
81e97e0d 553 //track - can be V0 according to TPC
afd25725 554 } else {
555 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
81e97e0d 556 delete t;
557 continue;
e50912db 558 }
afd25725 559 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
81e97e0d 560 delete t;
561 continue;
562 }
6c23ffed 563 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
81e97e0d 564 delete t;
565 continue;
566 }
e50912db 567 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
81e97e0d 568 delete t;
569 continue;
570 }
e43c066c 571 }
b9671574 572 t->SetReconstructed(kFALSE);
e43c066c 573 itsTracks.AddLast(t);
628e7bb0 574 fOriginal.AddLast(t);
e43c066c 575 }
576 } /* End Read ESD tracks */
577
578 itsTracks.Sort();
628e7bb0 579 fOriginal.Sort();
e43c066c 580 Int_t nentr=itsTracks.GetEntriesFast();
581 fTrackHypothesys.Expand(nentr);
81e97e0d 582 fBestHypothesys.Expand(nentr);
afd25725 583 MakeCoefficients(nentr);
513ab112 584 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
e43c066c 585 Int_t ntrk=0;
afd25725 586 // THE TWO TRACKING PASSES
e43c066c 587 for (fPass=0; fPass<2; fPass++) {
588 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
e50912db 589 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
e50912db 590 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
e43c066c 591 if (t==0) continue; //this track has been already tracked
1c97ce2f 592 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
b9671574 593 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
791f9a2a 594 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
afd25725 595 if (fConstraint[fPass]) {
596 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
597 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
598 }
e43c066c 599
600 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
4187a601 601 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
e43c066c 602 fI = 6;
603 ResetTrackToFollow(*t);
604 ResetBestTrack();
6518a6c5 605
e50912db 606 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
1c97ce2f 607
e43c066c 608
609 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
610 //
ae00569a 611 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
e43c066c 612 if (!besttrack) continue;
613 besttrack->SetLabel(tpcLabel);
614 // besttrack->CookdEdx();
615 CookdEdx(besttrack);
b9671574 616 besttrack->SetFakeRatio(1.);
e43c066c 617 CookLabel(besttrack,0.); //For comparison only
e43c066c 618 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
628e7bb0 619
628e7bb0 620 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
621
b9671574 622 t->SetReconstructed(kTRUE);
1c97ce2f 623 ntrk++;
4187a601 624 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
e43c066c 625 }
626 GetBestHypothesysMIP(itsTracks);
afd25725 627 } // end loop on the two tracking passes
e43c066c 628
cd90f0a2 629 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
630 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
628e7bb0 631 fAfterV0 = kTRUE;
628e7bb0 632 //
a0b53cc4 633 itsTracks.Clear();
e43c066c 634 //
635 Int_t entries = fTrackHypothesys.GetEntriesFast();
e50912db 636 for (Int_t ientry=0; ientry<entries; ientry++) {
e43c066c 637 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
638 if (array) array->Delete();
639 delete fTrackHypothesys.RemoveAt(ientry);
640 }
641
642 fTrackHypothesys.Delete();
4ffece29 643 entries = fBestHypothesys.GetEntriesFast();
644 for (Int_t ientry=0; ientry<entries; ientry++) {
645 TObjArray * array =(TObjArray*)fBestHypothesys.UncheckedAt(ientry);
646 if (array) array->Delete();
647 delete fBestHypothesys.RemoveAt(ientry);
648 }
81e97e0d 649 fBestHypothesys.Delete();
628e7bb0 650 fOriginal.Clear();
afd25725 651 delete [] fCoefficients;
652 fCoefficients=0;
e50912db 653 DeleteTrksMaterialLUT();
654
75a74b11 655 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
e50912db 656
657 fTrackingPhase="Default";
628e7bb0 658
e43c066c 659 return 0;
660}
44347160 661//------------------------------------------------------------------------
af885e0f 662Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
e43c066c 663 //--------------------------------------------------------------------
664 // This functions propagates reconstructed ITS tracks back
665 // The clusters must be loaded !
666 //--------------------------------------------------------------------
e50912db 667 fTrackingPhase="PropagateBack";
e43c066c 668 Int_t nentr=event->GetNumberOfTracks();
75a74b11 669 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
e43c066c 670
671 Int_t ntrk=0;
672 for (Int_t i=0; i<nentr; i++) {
673 AliESDtrack *esd=event->GetTrack(i);
674
675 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
676 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
677
763bdf3e 678 AliITStrackMI *t = new AliITStrackMI(*esd);
679
b9671574 680 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
e43c066c 681
682 ResetTrackToFollow(*t);
4fa7f7d1 683
67178f34 684 /*
9a5326b1 685 // propagate to vertex [SR, GSI 17.02.2003]
e43c066c 686 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
e50912db 687 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
688 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
689 fTrackToFollow.StartTimeIntegral();
690 // from vertex to outside pipe
691 CorrectForPipeMaterial(&fTrackToFollow,"outward");
67178f34 692 }*/
693 // Start time integral and add distance from current position to vertex
694 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
695 fTrackToFollow.GetXYZ(xyzTrk);
696 Double_t dst2 = 0.;
697 for (Int_t icoord=0; icoord<3; icoord++) {
698 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
699 dst2 += di*di;
e43c066c 700 }
67178f34 701 fTrackToFollow.StartTimeIntegral();
702 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
e43c066c 703
6c94f330 704 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
e50912db 705 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
706 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
e43c066c 707 delete t;
708 continue;
709 }
710 fTrackToFollow.SetLabel(t->GetLabel());
711 //fTrackToFollow.CookdEdx();
712 CookLabel(&fTrackToFollow,0.); //For comparison only
713 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
714 //UseClusters(&fTrackToFollow);
715 ntrk++;
716 }
717 delete t;
718 }
719
75a74b11 720 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
e43c066c 721
e50912db 722 fTrackingPhase="Default";
723
e43c066c 724 return 0;
725}
44347160 726//------------------------------------------------------------------------
af885e0f 727Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
e43c066c 728 //--------------------------------------------------------------------
729 // This functions refits ITS tracks using the
730 // "inward propagated" TPC tracks
731 // The clusters must be loaded !
732 //--------------------------------------------------------------------
e50912db 733 fTrackingPhase="RefitInward";
44396017 734
cd90f0a2 735 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
44396017 736
794d9013 737 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
738 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
739
e43c066c 740 Int_t nentr=event->GetNumberOfTracks();
75a74b11 741 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
e43c066c 742
743 Int_t ntrk=0;
744 for (Int_t i=0; i<nentr; i++) {
745 AliESDtrack *esd=event->GetTrack(i);
746
747 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
748 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
749 if (esd->GetStatus()&AliESDtrack::kTPCout)
750 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
751
763bdf3e 752 AliITStrackMI *t = new AliITStrackMI(*esd);
753
b9671574 754 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
e50912db 755 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
e43c066c 756 delete t;
757 continue;
758 }
759
760 ResetTrackToFollow(*t);
761 fTrackToFollow.ResetClusters();
762
ef3508c5 763 // ITS standalone tracks
764 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
6c94f330 765 fTrackToFollow.ResetCovariance(10.);
ef3508c5 766 // protection for loopers that can have parameters screwed up
767 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
768 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
769 delete t;
770 continue;
771 }
772 }
e43c066c 773
774 //Refitting...
58e8dc31 775 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
776 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
1a64bda6 777
4187a601 778 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
794d9013 779 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
4187a601 780 AliDebug(2," refit OK");
e43c066c 781 fTrackToFollow.SetLabel(t->GetLabel());
782 // fTrackToFollow.CookdEdx();
783 CookdEdx(&fTrackToFollow);
784
785 CookLabel(&fTrackToFollow,0.0); //For comparison only
786
9a5326b1 787 //The beam pipe
e50912db 788 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
ae00569a 789 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
9a5326b1 790 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
ae00569a 791 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
792 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
f7a1cc68 793 Double_t r[3]={0.,0.,0.};
9a5326b1 794 Double_t maxD=3.;
795 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
796 ntrk++;
e43c066c 797 }
798 }
799 delete t;
800 }
801
75a74b11 802 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
e43c066c 803
e50912db 804 fTrackingPhase="Default";
805
e43c066c 806 return 0;
807}
44347160 808//------------------------------------------------------------------------
e43c066c 809AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
810 //--------------------------------------------------------------------
811 // Return pointer to a given cluster
812 //--------------------------------------------------------------------
813 Int_t l=(index & 0xf0000000) >> 28;
814 Int_t c=(index & 0x0fffffff) >> 00;
815 return fgLayers[l].GetCluster(c);
816}
9a5326b1 817//------------------------------------------------------------------------
df29e9a4 818Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
9a5326b1 819 //--------------------------------------------------------------------
df29e9a4 820 // Get track space point with index i
9a5326b1 821 //--------------------------------------------------------------------
60066577 822
78b03929 823 Int_t l=(index & 0xf0000000) >> 28;
824 Int_t c=(index & 0x0fffffff) >> 00;
00a7cc50 825 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
df29e9a4 826 Int_t idet = cl->GetDetectorIndex();
4e65ff9e 827
60066577 828 Float_t xyz[3];
829 Float_t cov[6];
830 cl->GetGlobalXYZ(xyz);
831 cl->GetGlobalCov(cov);
832 p.SetXYZ(xyz, cov);
be77c8c9 833 p.SetCharge(cl->GetQ());
ebbcf1ea 834 p.SetDriftTime(cl->GetDriftTime());
b7bcc8ed 835 p.SetChargeRatio(cl->GetChargeRatio());
e3901fd8 836 p.SetClusterType(cl->GetClusterType());
ae079791 837 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
df29e9a4 838 switch (l) {
839 case 0:
ae079791 840 iLayer = AliGeomManager::kSPD1;
df29e9a4 841 break;
842 case 1:
ae079791 843 iLayer = AliGeomManager::kSPD2;
df29e9a4 844 break;
845 case 2:
ae079791 846 iLayer = AliGeomManager::kSDD1;
df29e9a4 847 break;
848 case 3:
ae079791 849 iLayer = AliGeomManager::kSDD2;
9a5326b1 850 break;
851 case 4:
852 iLayer = AliGeomManager::kSSD1;
853 break;
854 case 5:
855 iLayer = AliGeomManager::kSSD2;
856 break;
857 default:
858 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
859 break;
860 };
861 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
862 p.SetVolumeID((UShort_t)volid);
863 return kTRUE;
864}
865//------------------------------------------------------------------------
866Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
867 AliTrackPoint& p, const AliESDtrack *t) {
868 //--------------------------------------------------------------------
869 // Get track space point with index i
870 // (assign error estimated during the tracking)
871 //--------------------------------------------------------------------
872
873 Int_t l=(index & 0xf0000000) >> 28;
874 Int_t c=(index & 0x0fffffff) >> 00;
875 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
876 Int_t idet = cl->GetDetectorIndex();
4187a601 877
9a5326b1 878 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
879
880 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
881 Float_t detxy[2];
882 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
883 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
884 Double_t alpha = t->GetAlpha();
885 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
242bd643 886 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
9a5326b1 887 phi += alpha-det.GetPhi();
888 Float_t tgphi = TMath::Tan(phi);
889
890 Float_t tgl = t->GetTgl(); // tgl about const along track
891 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
892
d9ead1a0 893 Float_t errtrky,errtrkz,covyz;
8c139cf3 894 Bool_t addMisalErr=kFALSE;
d9ead1a0 895 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
9a5326b1 896
897 Float_t xyz[3];
898 Float_t cov[6];
899 cl->GetGlobalXYZ(xyz);
900 // cl->GetGlobalCov(cov);
901 Float_t pos[3] = {0.,0.,0.};
d9ead1a0 902 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
9a5326b1 903 tmpcl.GetGlobalCov(cov);
904
905 p.SetXYZ(xyz, cov);
be77c8c9 906 p.SetCharge(cl->GetQ());
ebbcf1ea 907 p.SetDriftTime(cl->GetDriftTime());
b7bcc8ed 908 p.SetChargeRatio(cl->GetChargeRatio());
e3901fd8 909 p.SetClusterType(cl->GetClusterType());
9a5326b1 910
911 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
912 switch (l) {
913 case 0:
914 iLayer = AliGeomManager::kSPD1;
915 break;
916 case 1:
917 iLayer = AliGeomManager::kSPD2;
918 break;
919 case 2:
920 iLayer = AliGeomManager::kSDD1;
921 break;
922 case 3:
923 iLayer = AliGeomManager::kSDD2;
df29e9a4 924 break;
925 case 4:
ae079791 926 iLayer = AliGeomManager::kSSD1;
df29e9a4 927 break;
928 case 5:
ae079791 929 iLayer = AliGeomManager::kSSD2;
df29e9a4 930 break;
931 default:
932 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
933 break;
934 };
ae079791 935 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
4187a601 936
df29e9a4 937 p.SetVolumeID((UShort_t)volid);
938 return kTRUE;
939}
44347160 940//------------------------------------------------------------------------
81e97e0d 941void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
e43c066c 942{
943 //--------------------------------------------------------------------
944 // Follow prolongation tree
945 //--------------------------------------------------------------------
81e97e0d 946 //
afd25725 947 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
948 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
949
2755f080 950
b9671574 951 AliESDtrack * esd = otrack->GetESDtrack();
afd25725 952 if (esd->GetV0Index(0)>0) {
81e97e0d 953 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
afd25725 954 // mapping of ESD track is different as ITS track in Containers
81e97e0d 955 // Need something more stable
44347160 956 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
81e97e0d 957 for (Int_t i=0;i<3;i++){
958 Int_t index = esd->GetV0Index(i);
959 if (index==0) break;
d6a49f20 960 AliESDv0 * vertex = fEsd->GetV0(index);
81e97e0d 961 if (vertex->GetStatus()<0) continue; // rejected V0
962 //
afd25725 963 if (esd->GetSign()>0) {
964 vertex->SetIndex(0,esdindex);
965 } else {
966 vertex->SetIndex(1,esdindex);
967 }
81e97e0d 968 }
969 }
970 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
971 if (!bestarray){
972 bestarray = new TObjArray(5);
4ffece29 973 bestarray->SetOwner();
81e97e0d 974 fBestHypothesys.AddAt(bestarray,esdindex);
975 }
e43c066c 976
81e97e0d 977 //
e43c066c 978 //setup tree of the prolongations
979 //
15dd636f 980 static AliITStrackMI tracks[7][100];
981 AliITStrackMI *currenttrack;
982 static AliITStrackMI currenttrack1;
983 static AliITStrackMI currenttrack2;
984 static AliITStrackMI backuptrack;
e43c066c 985 Int_t ntracks[7];
986 Int_t nindexes[7][100];
987 Float_t normalizedchi2[100];
988 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
b9671574 989 otrack->SetNSkipped(0);
15dd636f 990 new (&(tracks[6][0])) AliITStrackMI(*otrack);
e43c066c 991 ntracks[6]=1;
81e97e0d 992 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
ae00569a 993 Int_t modstatus = 1; // found
994 Float_t xloc,zloc;
e43c066c 995 //
996 //
997 // follow prolongations
e50912db 998 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
4187a601 999 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
2755f080 1000 fI = ilayer;
e43c066c 1001 //
2755f080 1002 AliITSlayer &layer=fgLayers[ilayer];
1003 Double_t r = layer.GetR();
e43c066c 1004 ntracks[ilayer]=0;
1005 //
1006 //
2755f080 1007 Int_t nskipped=0;
e43c066c 1008 Float_t nused =0;
e50912db 1009 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
e43c066c 1010 //set current track
1011 if (ntracks[ilayer]>=100) break;
b9671574 1012 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1013 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
e43c066c 1014 if (ntracks[ilayer]>15+ilayer){
b9671574 1015 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1016 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
e43c066c 1017 }
1018
15dd636f 1019 new(&currenttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
44347160 1020
afd25725 1021 // material between SSD and SDD, SDD and SPD
e50912db 1022 if (ilayer==3)
1023 if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
1024 if (ilayer==1)
1025 if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
8602c008 1026
2755f080 1027 // detector number
8602c008 1028 Double_t phi,z;
1029 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
e43c066c 1030 Int_t idet=layer.FindDetectorIndex(phi,z);
e50912db 1031
1032 Double_t trackGlobXYZ1[3];
d1181f10 1033 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
e50912db 1034
2755f080 1035 // Get the budget to the primary vertex for the current track being prolonged
1036 Double_t budgetToPrimVertex = GetEffectiveThickness();
1037
1038 // check if we allow a prolongation without point
6518a6c5 1039 Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
2755f080 1040 if (skip) {
1041 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
6518a6c5 1042 // propagate to the layer radius
d1181f10 1043 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
9be1d1c7 1044 if(!vtrack->Propagate(xToGo)) continue;
2755f080 1045 // apply correction for material of the current layer
1046 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1047 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
f8720bda 1048 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
9fd412a8 1049 vtrack->SetClIndex(ilayer,-1);
ae00569a 1050 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
d1181f10 1051 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1052 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1053 }
333d86cb 1054 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
2755f080 1055 ntracks[ilayer]++;
6518a6c5 1056 continue;
2755f080 1057 }
1058
1059 // track outside layer acceptance in z
1060 if (idet<0) continue;
1061
1062 //propagate to the intersection with the detector plane
1063 const AliITSdetector &det=layer.GetDetector(idet);
1064 new(&currenttrack2) AliITStrackMI(currenttrack1);
44347160 1065 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
d1181f10 1066 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
e43c066c 1067 currenttrack1.SetDetectorIndex(idet);
1068 currenttrack2.SetDetectorIndex(idet);
d1181f10 1069 if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
6518a6c5 1070
afd25725 1071 //***************
1c97ce2f 1072 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
e43c066c 1073 //
2755f080 1074 // road in global (rphi,z) [i.e. in tracking ref. system]
1c97ce2f 1075 Double_t zmin,zmax,ymin,ymax;
1076 if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
2755f080 1077
afd25725 1078 // select clusters in road
e43c066c 1079 layer.SelectClusters(zmin,zmax,ymin,ymax);
afd25725 1080 //********************
44347160 1081
afd25725 1082 // Define criteria for track-cluster association
44347160 1083 Double_t msz = currenttrack1.GetSigmaZ2() +
1084 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1085 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1086 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1087 Double_t msy = currenttrack1.GetSigmaY2() +
1088 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1089 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1090 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1091 if (constrain) {
1092 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1093 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1094 } else {
1095 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1096 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1097 }
1098 msz = 1./msz; // 1/RoadZ^2
1099 msy = 1./msy; // 1/RoadY^2
1c97ce2f 1100
e43c066c 1101 //
e43c066c 1102 //
afd25725 1103 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
e43c066c 1104 //
afd25725 1105 const AliITSRecPoint *cl=0;
1106 Int_t clidx=-1;
1107 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
6518a6c5 1108 Bool_t deadzoneSPD=kFALSE;
e43c066c 1109 currenttrack = &currenttrack1;
6518a6c5 1110
ae00569a 1111 // check if the road contains a dead zone
23197852 1112 Bool_t noClusters = kFALSE;
1113 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
4187a601 1114 if (noClusters) AliDebug(2,"no clusters in road");
1c97ce2f 1115 Double_t dz=0.5*(zmax-zmin);
1116 Double_t dy=0.5*(ymax-ymin);
1117 Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters);
4187a601 1118 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
ae00569a 1119 // create a prolongation without clusters (check also if there are no clusters in the road)
1120 if (dead ||
23197852 1121 (noClusters &&
ae00569a 1122 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
6518a6c5 1123 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
9fd412a8 1124 updatetrack->SetClIndex(ilayer,-1);
ae00569a 1125 if (dead==0) {
1126 modstatus = 5; // no cls in road
1127 } else if (dead==1) {
1128 modstatus = 7; // holes in z in SPD
f8720bda 1129 } else if (dead==2 || dead==3 || dead==4) {
ae00569a 1130 modstatus = 2; // dead from OCDB
1131 }
1132 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
6518a6c5 1133 // apply correction for material of the current layer
1134 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1135 if (constrain) { // apply vertex constrain
1136 updatetrack->SetConstrain(constrain);
1137 Bool_t isPrim = kTRUE;
1138 if (ilayer<4) { // check that it's close to the vertex
1139 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1140 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1141 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1142 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1143 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1144 }
333d86cb 1145 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
6518a6c5 1146 }
f8720bda 1147 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
ae00569a 1148 if (dead) {
ae00569a 1149 if (dead==1) { // dead zone at z=0,+-7cm in SPD
f8720bda 1150 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
ae00569a 1151 deadzoneSPD=kTRUE;
f8720bda 1152 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1153 updatetrack->SetDeadZoneProbability(ilayer,1.);
1154 } else if (dead==4) { // at least a single dead channel from OCDB
1155 updatetrack->SetDeadZoneProbability(ilayer,0.);
1156 }
6518a6c5 1157 }
1158 ntracks[ilayer]++;
1159 }
1160
ae00569a 1161 clidx=-1;
6518a6c5 1162 // loop over clusters in the road
afd25725 1163 while ((cl=layer.GetNextCluster(clidx))!=0) {
e43c066c 1164 if (ntracks[ilayer]>95) break; //space for skipped clusters
afd25725 1165 Bool_t changedet =kFALSE;
88e3c173 1166 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
a5d0fae6 1167 Int_t idetc=cl->GetDetectorIndex();
afd25725 1168
a5d0fae6 1169 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1c97ce2f 1170 // take into account misalignment (bring track to real detector plane)
1171 Double_t xTrOrig = currenttrack->GetX();
9be1d1c7 1172 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
afd25725 1173 // a first cut on track-cluster distance
1174 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1175 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1c97ce2f 1176 { // cluster not associated to track
4187a601 1177 AliDebug(2,"not associated");
421c4ba7 1178 // MvL: added here as well
1179 // bring track back to ideal detector plane
1180 currenttrack->Propagate(xTrOrig);
1c97ce2f 1181 continue;
1182 }
1183 // bring track back to ideal detector plane
9be1d1c7 1184 if (!currenttrack->Propagate(xTrOrig)) continue;
afd25725 1185 } else { // have to move track to cluster's detector
a5d0fae6 1186 const AliITSdetector &detc=layer.GetDetector(idetc);
afd25725 1187 // a first cut on track-cluster distance
a5d0fae6 1188 Double_t y;
1c97ce2f 1189 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
afd25725 1190 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1191 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
44347160 1192 continue; // cluster not associated to track
e43c066c 1193 //
15dd636f 1194 new (&backuptrack) AliITStrackMI(currenttrack2);
afd25725 1195 changedet = kTRUE;
e43c066c 1196 currenttrack =&currenttrack2;
a5d0fae6 1197 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
15dd636f 1198 new (currenttrack) AliITStrackMI(backuptrack);
afd25725 1199 changedet = kFALSE;
e43c066c 1200 continue;
1201 }
a5d0fae6 1202 currenttrack->SetDetectorIndex(idetc);
e50912db 1203 // Get again the budget to the primary vertex
1204 // for the current track being prolonged, if had to change detector
1205 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
e43c066c 1206 }
1207
afd25725 1208 // calculate track-clusters chi2
1209 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1210 // chi2 cut
4187a601 1211 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
afd25725 1212 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
88e3c173 1213 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
e43c066c 1214 if (ntracks[ilayer]>=100) continue;
15dd636f 1215 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
9fd412a8 1216 updatetrack->SetClIndex(ilayer,-1);
afd25725 1217 if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
ae00569a 1218
88e3c173 1219 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1c97ce2f 1220 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
4187a601 1221 AliDebug(2,"update failed");
1c97ce2f 1222 continue;
1223 }
09cf9d66 1224 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
ae00569a 1225 modstatus = 1; // found
6518a6c5 1226 } else { // virtual cluster in dead zone
b9671574 1227 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
f8720bda 1228 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
ae00569a 1229 modstatus = 7; // holes in z in SPD
1230 }
1231
1232 if (changedet) {
1233 Float_t xlocnewdet,zlocnewdet;
d1181f10 1234 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1235 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1236 }
ae00569a 1237 } else {
1238 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
e43c066c 1239 }
afd25725 1240 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1241
1242 // apply correction for material of the current layer
e50912db 1243 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1244
afd25725 1245 if (constrain) { // apply vertex constrain
b9671574 1246 updatetrack->SetConstrain(constrain);
e43c066c 1247 Bool_t isPrim = kTRUE;
afd25725 1248 if (ilayer<4) { // check that it's close to the vertex
b9671574 1249 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
afd25725 1250 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1251 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1252 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1253 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
e43c066c 1254 }
333d86cb 1255 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
e43c066c 1256 } //apply vertex constrain
1257 ntracks[ilayer]++;
afd25725 1258 } // create new hypothesis
4187a601 1259 else {
1260 AliDebug(2,"chi2 too large");
1261 }
1262
2755f080 1263 } // loop over possible prolongations
1264
1265 // allow one prolongation without clusters
88e3c173 1266 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
15dd636f 1267 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
2755f080 1268 // apply correction for material of the current layer
1269 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
9fd412a8 1270 vtrack->SetClIndex(ilayer,-1);
ae00569a 1271 modstatus = 3; // skipped
1272 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
333d86cb 1273 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
b9671574 1274 vtrack->IncrementNSkipped();
e43c066c 1275 ntracks[ilayer]++;
1276 }
1277
6262dd3d 1278
e50912db 1279 } // loop over tracks in layer ilayer+1
afd25725 1280
1281 //loop over track candidates for the current layer
e43c066c 1282 //
1283 //
1284 Int_t accepted=0;
1285
afd25725 1286 Int_t golden=0;
e43c066c 1287 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1288 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
afd25725 1289 if (normalizedchi2[itrack] <
1290 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
44347160 1291 if (ilayer>4) {
1292 accepted++;
1293 } else {
1294 if (constrain) { // constrain
afd25725 1295 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1296 accepted++;
44347160 1297 } else { // !constrain
afd25725 1298 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1299 accepted++;
44347160 1300 }
e43c066c 1301 }
1302 }
afd25725 1303 // sort tracks by increasing normalized chi2
1304 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
e43c066c 1305 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
afd25725 1306 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
e43c066c 1307 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
afd25725 1308 } // end loop over layers
1309
afd25725 1310
1311 //
1312 // Now select tracks to be kept
1313 //
44347160 1314 Int_t max = constrain ? 20 : 5;
e43c066c 1315
afd25725 1316 // tracks that reach layer 0 (SPD inner)
44347160 1317 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
15dd636f 1318 AliITStrackMI & track= tracks[0][nindexes[0][i]];
628e7bb0 1319 if (track.GetNumberOfClusters()<2) continue;
afd25725 1320 if (!constrain && track.GetNormChi2(0) >
6518a6c5 1321 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1322 continue;
1323 }
15dd636f 1324 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
e43c066c 1325 }
afd25725 1326
1327 // tracks that reach layer 1 (SPD outer)
628e7bb0 1328 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
15dd636f 1329 AliITStrackMI & track= tracks[1][nindexes[1][i]];
e43c066c 1330 if (track.GetNumberOfClusters()<4) continue;
afd25725 1331 if (!constrain && track.GetNormChi2(1) >
1332 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
b9671574 1333 if (constrain) track.IncrementNSkipped();
81e97e0d 1334 if (!constrain) {
b9671574 1335 track.SetD(0,track.GetD(GetX(),GetY()));
afd25725 1336 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
b9671574 1337 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1338 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
e43c066c 1339 }
1340 }
15dd636f 1341 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
e43c066c 1342 }
afd25725 1343
6518a6c5 1344 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
81e97e0d 1345 if (!constrain){
628e7bb0 1346 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
15dd636f 1347 AliITStrackMI & track= tracks[2][nindexes[2][i]];
628e7bb0 1348 if (track.GetNumberOfClusters()<3) continue;
9cd19c63 1349 if (track.GetNormChi2(2) >
afd25725 1350 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
9cd19c63 1351 track.SetD(0,track.GetD(GetX(),GetY()));
1352 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1353 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1354 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
e43c066c 1355 }
15dd636f 1356 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
e43c066c 1357 }
1358 }
81e97e0d 1359
afd25725 1360 if (!constrain) {
81e97e0d 1361 //
afd25725 1362 // register best track of each layer - important for V0 finder
81e97e0d 1363 //
1364 for (Int_t ilayer=0;ilayer<5;ilayer++){
1365 if (ntracks[ilayer]==0) continue;
1366 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1367 if (track.GetNumberOfClusters()<1) continue;
1368 CookLabel(&track,0);
1369 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1370 }
1371 }
1372 //
1373 // update TPC V0 information
1374 //
b9671574 1375 if (otrack->GetESDtrack()->GetV0Index(0)>0){
81e97e0d 1376 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1377 for (Int_t i=0;i<3;i++){
b9671574 1378 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
81e97e0d 1379 if (index==0) break;
44347160 1380 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
81e97e0d 1381 if (vertex->GetStatus()<0) continue; // rejected V0
1382 //
6c94f330 1383 if (otrack->GetSign()>0) {
81e97e0d 1384 vertex->SetIndex(0,esdindex);
1385 }
1386 else{
1387 vertex->SetIndex(1,esdindex);
1388 }
1389 //find nearest layer with track info
b75d63a7 1390 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1391 Int_t nearestold = GetNearestLayer(xrp); //I.B.
81e97e0d 1392 Int_t nearest = nearestold;
9cd19c63 1393 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
81e97e0d 1394 if (ntracks[nearest]==0){
1395 nearest = ilayer;
1396 }
1397 }
1398 //
1399 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1400 if (nearestold<5&&nearest<5){
b9671574 1401 Bool_t accept = track.GetNormChi2(nearest)<10;
81e97e0d 1402 if (accept){
6c94f330 1403 if (track.GetSign()>0) {
b75d63a7 1404 vertex->SetParamP(track);
81e97e0d 1405 vertex->Update(fprimvertex);
44347160 1406 //vertex->SetIndex(0,track.fESDtrack->GetID());
81e97e0d 1407 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1408 }else{
b75d63a7 1409 vertex->SetParamN(track);
81e97e0d 1410 vertex->Update(fprimvertex);
1411 //vertex->SetIndex(1,track.fESDtrack->GetID());
1412 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1413 }
1414 vertex->SetStatus(vertex->GetStatus()+1);
1415 }else{
44347160 1416 //vertex->SetStatus(-2); // reject V0 - not enough clusters
81e97e0d 1417 }
1418 }
81e97e0d 1419 }
afd25725 1420 }
1421
e43c066c 1422}
44347160 1423//------------------------------------------------------------------------
e43c066c 1424AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1425{
1426 //--------------------------------------------------------------------
1427 //
1428 //
1429 return fgLayers[layer];
1430}
44347160 1431//------------------------------------------------------------------------
8221b41b 1432AliITStrackerMI::AliITSlayer::AliITSlayer():
1433fR(0),
1434fPhiOffset(0),
1435fNladders(0),
1436fZOffset(0),
1437fNdetectors(0),
1438fDetectors(0),
1439fN(0),
1440fDy5(0),
1441fDy10(0),
1442fDy20(0),
1443fClustersCs(0),
1444fClusterIndexCs(0),
1445fYcs(0),
1446fZcs(0),
1447fNcs(0),
1448fCurrentSlice(-1),
d9ead1a0 1449fZmin(0),
8221b41b 1450fZmax(0),
1451fYmin(0),
1452fYmax(0),
1453fI(0),
1454fImax(0),
1455fSkip(0),
1456fAccepted(0),
d9ead1a0 1457fRoad(0),
1458fMaxSigmaClY(0),
1459fMaxSigmaClZ(0),
1460fNMaxSigmaCl(3)
1461{
e43c066c 1462 //--------------------------------------------------------------------
1463 //default AliITSlayer constructor
1464 //--------------------------------------------------------------------
e50912db 1465 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
e43c066c 1466 fClusterWeight[i]=0;
1467 fClusterTracks[0][i]=-1;
1468 fClusterTracks[1][i]=-1;
1469 fClusterTracks[2][i]=-1;
9cd19c63 1470 fClusterTracks[3][i]=-1;
1471 fY[i]=0;
1472 fZ[i]=0;
1473 }
1474 fYB[0]=0;
1475 fYB[1]=0;
1476
1477 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
1478 for (Int_t j1=0; j1<6; j1++) {
1479 fClusters5[j1][j]=0;
1480 fClusterIndex5[j1][j]=-1;
1481 fY5[j1][j]=0;
1482 fZ5[j1][j]=0;
1483 fN5[j1]=0;
1484 fBy5[j1][0]=0;
1485 fBy5[j1][1]=0;
1486 }
1487 }
1488
1489 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
1490 for (Int_t j1=0; j1<11; j1++) {
1491 fClusters10[j1][j]=0;
1492 fClusterIndex10[j1][j]=-1;
1493 fY10[j1][j]=0;
1494 fZ10[j1][j]=0;
1495 fN10[j1]=0;
1496 fBy10[j1][0]=0;
1497 fBy10[j1][1]=0;
1498 }
e43c066c 1499 }
9cd19c63 1500
1501 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
1502 for (Int_t j1=0; j1<21; j1++) {
1503 fClusters20[j1][j]=0;
1504 fClusterIndex20[j1][j]=-1;
1505 fY20[j1][j]=0;
1506 fZ20[j1][j]=0;
1507 fN20[j1]=0;
1508 fBy20[j1][0]=0;
1509 fBy20[j1][1]=0;
1510 }
1511 }
e8c2bf0d 1512 for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
1513 fClusters[i]=NULL;
1514 fClusterIndex[i]=0;
1515 }
e43c066c 1516}
44347160 1517//------------------------------------------------------------------------
e43c066c 1518AliITStrackerMI::AliITSlayer::
8221b41b 1519AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1520fR(r),
1521fPhiOffset(p),
1522fNladders(nl),
1523fZOffset(z),
1524fNdetectors(nd),
1525fDetectors(0),
1526fN(0),
1527fDy5(0),
1528fDy10(0),
1529fDy20(0),
1530fClustersCs(0),
1531fClusterIndexCs(0),
1532fYcs(0),
1533fZcs(0),
1534fNcs(0),
1535fCurrentSlice(-1),
d9ead1a0 1536fZmin(0),
8221b41b 1537fZmax(0),
1538fYmin(0),
1539fYmax(0),
1540fI(0),
1541fImax(0),
1542fSkip(0),
1543fAccepted(0),
d9ead1a0 1544fRoad(0),
1545fMaxSigmaClY(0),
1546fMaxSigmaClZ(0),
1547fNMaxSigmaCl(3) {
e43c066c 1548 //--------------------------------------------------------------------
1549 //main AliITSlayer constructor
1550 //--------------------------------------------------------------------
e43c066c 1551 fDetectors=new AliITSdetector[fNladders*fNdetectors];
afd25725 1552 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
9cd19c63 1553
1554 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1555 fClusterWeight[i]=0;
1556 fClusterTracks[0][i]=-1;
1557 fClusterTracks[1][i]=-1;
1558 fClusterTracks[2][i]=-1;
1559 fClusterTracks[3][i]=-1;
1560 fY[i]=0;
1561 fZ[i]=0;
1562 }
1563
1564 fYB[0]=0;
1565 fYB[1]=0;
1566
1567 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
1568 for (Int_t j1=0; j1<6; j1++) {
1569 fClusters5[j1][j]=0;
1570 fClusterIndex5[j1][j]=-1;
1571 fY5[j1][j]=0;
1572 fZ5[j1][j]=0;
1573 fN5[j1]=0;
1574 fBy5[j1][0]=0;
1575 fBy5[j1][1]=0;
1576 }
1577 }
1578
1579 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
1580 for (Int_t j1=0; j1<11; j1++) {
1581 fClusters10[j1][j]=0;
1582 fClusterIndex10[j1][j]=-1;
1583 fY10[j1][j]=0;
1584 fZ10[j1][j]=0;
1585 fN10[j1]=0;
1586 fBy10[j1][0]=0;
1587 fBy10[j1][1]=0;
1588 }
1589 }
1590
1591 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
1592 for (Int_t j1=0; j1<21; j1++) {
1593 fClusters20[j1][j]=0;
1594 fClusterIndex20[j1][j]=-1;
1595 fY20[j1][j]=0;
1596 fZ20[j1][j]=0;
1597 fN20[j1]=0;
1598 fBy20[j1][0]=0;
1599 fBy20[j1][1]=0;
1600 }
1601 }
e8c2bf0d 1602 for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
1603 fClusters[i]=NULL;
1604 fClusterIndex[i]=0;
1605 }
e43c066c 1606}
225cff59 1607/*
44347160 1608//------------------------------------------------------------------------
8221b41b 1609AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1610fR(layer.fR),
1611fPhiOffset(layer.fPhiOffset),
1612fNladders(layer.fNladders),
1613fZOffset(layer.fZOffset),
1614fNdetectors(layer.fNdetectors),
1615fDetectors(layer.fDetectors),
1616fN(layer.fN),
1617fDy5(layer.fDy5),
1618fDy10(layer.fDy10),
1619fDy20(layer.fDy20),
1620fClustersCs(layer.fClustersCs),
1621fClusterIndexCs(layer.fClusterIndexCs),
1622fYcs(layer.fYcs),
1623fZcs(layer.fZcs),
1624fNcs(layer.fNcs),
1625fCurrentSlice(layer.fCurrentSlice),
d9ead1a0 1626fZmin(layer.fZmin),
8221b41b 1627fZmax(layer.fZmax),
1628fYmin(layer.fYmin),
1629fYmax(layer.fYmax),
1630fI(layer.fI),
1631fImax(layer.fImax),
1632fSkip(layer.fSkip),
1633fAccepted(layer.fAccepted),
d9ead1a0 1634fRoad(layer.fRoad),
1635fMaxSigmaClY(layer.fMaxSigmaClY),
1636fMaxSigmaClZ(layer.fMaxSigmaClZ),
1637fNMaxSigmaCl(layer.fNMaxSigmaCl)
1638{
8221b41b 1639 //Copy constructor
1640}
225cff59 1641*/
44347160 1642//------------------------------------------------------------------------
e43c066c 1643AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1644 //--------------------------------------------------------------------
1645 // AliITSlayer destructor
1646 //--------------------------------------------------------------------
23197852 1647 delete [] fDetectors;
e43c066c 1648 for (Int_t i=0; i<fN; i++) delete fClusters[i];
e50912db 1649 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
e43c066c 1650 fClusterWeight[i]=0;
1651 fClusterTracks[0][i]=-1;
1652 fClusterTracks[1][i]=-1;
1653 fClusterTracks[2][i]=-1;
1654 fClusterTracks[3][i]=-1;
1655 }
1656}
44347160 1657//------------------------------------------------------------------------
e43c066c 1658void AliITStrackerMI::AliITSlayer::ResetClusters() {
1659 //--------------------------------------------------------------------
1660 // This function removes loaded clusters
1661 //--------------------------------------------------------------------
1662 for (Int_t i=0; i<fN; i++) delete fClusters[i];
e50912db 1663 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
e43c066c 1664 fClusterWeight[i]=0;
1665 fClusterTracks[0][i]=-1;
1666 fClusterTracks[1][i]=-1;
1667 fClusterTracks[2][i]=-1;
1668 fClusterTracks[3][i]=-1;
1669 }
1670
1671 fN=0;
1672 fI=0;
1673}
44347160 1674//------------------------------------------------------------------------
e43c066c 1675void AliITStrackerMI::AliITSlayer::ResetWeights() {
1676 //--------------------------------------------------------------------
1677 // This function reset weights of the clusters
1678 //--------------------------------------------------------------------
e50912db 1679 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
e43c066c 1680 fClusterWeight[i]=0;
1681 fClusterTracks[0][i]=-1;
1682 fClusterTracks[1][i]=-1;
1683 fClusterTracks[2][i]=-1;
1684 fClusterTracks[3][i]=-1;
1685 }
1686 for (Int_t i=0; i<fN;i++) {
00a7cc50 1687 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
e43c066c 1688 if (cl&&cl->IsUsed()) cl->Use();
1689 }
1690
1691}
44347160 1692//------------------------------------------------------------------------
e43c066c 1693void AliITStrackerMI::AliITSlayer::ResetRoad() {
1694 //--------------------------------------------------------------------
1695 // This function calculates the road defined by the cluster density
1696 //--------------------------------------------------------------------
1697 Int_t n=0;
1698 for (Int_t i=0; i<fN; i++) {
1699 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1700 }
afd25725 1701 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
e43c066c 1702}
44347160 1703//------------------------------------------------------------------------
afd25725 1704Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
e43c066c 1705 //--------------------------------------------------------------------
1706 //This function adds a cluster to this layer
1707 //--------------------------------------------------------------------
e50912db 1708 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
e43c066c 1709 return 1;
1710 }
1711 fCurrentSlice=-1;
afd25725 1712 fClusters[fN]=cl;
1d4090b7 1713 fN++;
afd25725 1714 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
d9ead1a0 1715 //AD
1716 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1717 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1718 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1719 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1720 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1721 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1722 //AD
1723 /*
afd25725 1724 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1725 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1726 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1727 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
d9ead1a0 1728 */
e43c066c 1729 return 0;
1730}
44347160 1731//------------------------------------------------------------------------
e43c066c 1732void AliITStrackerMI::AliITSlayer::SortClusters()
1733{
1734 //
1735 //sort clusters
1736 //
00a7cc50 1737 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1d4090b7 1738 Float_t *z = new Float_t[fN];
1739 Int_t * index = new Int_t[fN];
1740 //
d9ead1a0 1741 fMaxSigmaClY=0.; //AD
1742 fMaxSigmaClZ=0.; //AD
1743
1d4090b7 1744 for (Int_t i=0;i<fN;i++){
1745 z[i] = fClusters[i]->GetZ();
d9ead1a0 1746 // save largest errors in y and z for this layer
1747 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1748 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1d4090b7 1749 }
1750 TMath::Sort(fN,z,index,kFALSE);
1751 for (Int_t i=0;i<fN;i++){
1752 clusters[i] = fClusters[index[i]];
1753 }
1754 //
1755 for (Int_t i=0;i<fN;i++){
1756 fClusters[i] = clusters[i];
1757 fZ[i] = fClusters[i]->GetZ();
1758 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1759 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1760 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1761 fY[i] = y;
1762 }
1763 delete[] index;
1764 delete[] z;
1765 delete[] clusters;
1766 //
1767
e43c066c 1768 fYB[0]=10000000;
1769 fYB[1]=-10000000;
1770 for (Int_t i=0;i<fN;i++){
1771 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1772 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1773 fClusterIndex[i] = i;
1774 }
1775 //
1776 // fill slices
1777 fDy5 = (fYB[1]-fYB[0])/5.;
1778 fDy10 = (fYB[1]-fYB[0])/10.;
1779 fDy20 = (fYB[1]-fYB[0])/20.;
1780 for (Int_t i=0;i<6;i++) fN5[i] =0;
1781 for (Int_t i=0;i<11;i++) fN10[i]=0;
1782 for (Int_t i=0;i<21;i++) fN20[i]=0;
1783 //
1784 for (Int_t i=0;i<6;i++) {fBy5[i][0] = fYB[0]+(i-0.75)*fDy5; fBy5[i][1] = fYB[0]+(i+0.75)*fDy5;}
1785 for (Int_t i=0;i<11;i++) {fBy10[i][0] = fYB[0]+(i-0.75)*fDy10; fBy10[i][1] = fYB[0]+(i+0.75)*fDy10;}
1786 for (Int_t i=0;i<21;i++) {fBy20[i][0] = fYB[0]+(i-0.75)*fDy20; fBy20[i][1] = fYB[0]+(i+0.75)*fDy20;}
1787 //
1788 //
1d4090b7 1789 for (Int_t i=0;i<fN;i++)
1790 for (Int_t irot=-1;irot<=1;irot++){
1791 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1792 // slice 5
1793 for (Int_t slice=0; slice<6;slice++){
e50912db 1794 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1d4090b7 1795 fClusters5[slice][fN5[slice]] = fClusters[i];
1796 fY5[slice][fN5[slice]] = curY;
1797 fZ5[slice][fN5[slice]] = fZ[i];
1798 fClusterIndex5[slice][fN5[slice]]=i;
1799 fN5[slice]++;
1800 }
e43c066c 1801 }
1d4090b7 1802 // slice 10
1803 for (Int_t slice=0; slice<11;slice++){
e50912db 1804 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1d4090b7 1805 fClusters10[slice][fN10[slice]] = fClusters[i];
1806 fY10[slice][fN10[slice]] = curY;
1807 fZ10[slice][fN10[slice]] = fZ[i];
1808 fClusterIndex10[slice][fN10[slice]]=i;
1809 fN10[slice]++;
1810 }
e43c066c 1811 }
1d4090b7 1812 // slice 20
1813 for (Int_t slice=0; slice<21;slice++){
e50912db 1814 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1d4090b7 1815 fClusters20[slice][fN20[slice]] = fClusters[i];
1816 fY20[slice][fN20[slice]] = curY;
1817 fZ20[slice][fN20[slice]] = fZ[i];
1818 fClusterIndex20[slice][fN20[slice]]=i;
1819 fN20[slice]++;
1820 }
1821 }
e43c066c 1822 }
1d4090b7 1823
1824 //
1825 // consistency check
1826 //
1827 for (Int_t i=0;i<fN-1;i++){
1828 if (fZ[i]>fZ[i+1]){
ae00569a 1829 printf("Bug\n");
e43c066c 1830 }
1831 }
1d4090b7 1832 //
1833 for (Int_t slice=0;slice<21;slice++)
1834 for (Int_t i=0;i<fN20[slice]-1;i++){
1835 if (fZ20[slice][i]>fZ20[slice][i+1]){
ae00569a 1836 printf("Bug\n");
1d4090b7 1837 }
1838 }
1839
1840
e43c066c 1841}
44347160 1842//------------------------------------------------------------------------
e43c066c 1843Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1844 //--------------------------------------------------------------------
1845 // This function returns the index of the nearest cluster
1846 //--------------------------------------------------------------------
1847 Int_t ncl=0;
1848 const Float_t *zcl;
1849 if (fCurrentSlice<0) {
1850 ncl = fN;
1851 zcl = fZ;
1852 }
1853 else{
1854 ncl = fNcs;
1855 zcl = fZcs;;
1856 }
1857
1858 if (ncl==0) return 0;
1859 Int_t b=0, e=ncl-1, m=(b+e)/2;
1860 for (; b<e; m=(b+e)/2) {
1861 // if (z > fClusters[m]->GetZ()) b=m+1;
1862 if (z > zcl[m]) b=m+1;
1863 else e=m;
1864 }
1865 return m;
1866}
44347160 1867//------------------------------------------------------------------------
1c97ce2f 1868Bool_t AliITStrackerMI::ComputeRoad(AliITStrackMI* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1869 //--------------------------------------------------------------------
1870 // This function computes the rectangular road for this track
1871 //--------------------------------------------------------------------
1872
1873
1874 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1875 // take into account the misalignment: propagate track to misaligned detector plane
1876 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1877
1878 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1879 TMath::Sqrt(track->GetSigmaZ2() +
1880 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1881 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1882 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1883 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1884 TMath::Sqrt(track->GetSigmaY2() +
1885 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1886 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1887 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1888
1889 // track at boundary between detectors, enlarge road
1890 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1891 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1892 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1893 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1894 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1895 Float_t tgl = TMath::Abs(track->GetTgl());
1896 if (tgl > 1.) tgl=1.;
1897 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1898 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1899 Float_t snp = TMath::Abs(track->GetSnp());
1900 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1901 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1902 } // boundary
1903
1904 // add to the road a term (up to 2-3 mm) to deal with misalignments
1905 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1906 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1907
1908 Double_t r = fgLayers[ilayer].GetR();
1909 zmin = track->GetZ() - dz;
1910 zmax = track->GetZ() + dz;
1911 ymin = track->GetY() + r*det.GetPhi() - dy;
1912 ymax = track->GetY() + r*det.GetPhi() + dy;
1913
1914 // bring track back to idead detector plane
1915 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1916
1917 return kTRUE;
1918}
1919//------------------------------------------------------------------------
e43c066c 1920void AliITStrackerMI::AliITSlayer::
1921SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1922 //--------------------------------------------------------------------
1923 // This function sets the "window"
1924 //--------------------------------------------------------------------
1925
628e7bb0 1926 Double_t circle=2*TMath::Pi()*fR;
d9ead1a0 1927 fYmin = ymin;
1928 fYmax = ymax;
1929 fZmin = zmin;
1930 fZmax = zmax;
1931 // AD
1932 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1933 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1934 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1935 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1936 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1937
e43c066c 1938 Float_t ymiddle = (fYmin+fYmax)*0.5;
ae00569a 1939 if (ymiddle<fYB[0]) {
1940 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1941 } else if (ymiddle>fYB[1]) {
1942 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
e43c066c 1943 }
ae00569a 1944
e43c066c 1945 //
1946 fCurrentSlice =-1;
1947 // defualt take all
1948 fClustersCs = fClusters;
1949 fClusterIndexCs = fClusterIndex;
1950 fYcs = fY;
1951 fZcs = fZ;
1952 fNcs = fN;
628e7bb0 1953 //
e43c066c 1954 //is in 20 slice?
1955 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1956 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1957 if (slice<0) slice=0;
1958 if (slice>20) slice=20;
1959 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1960 if (isOK) {
1961 fCurrentSlice=slice;
1962 fClustersCs = fClusters20[fCurrentSlice];
1963 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1964 fYcs = fY20[fCurrentSlice];
1965 fZcs = fZ20[fCurrentSlice];
1966 fNcs = fN20[fCurrentSlice];
1967 }
1968 }
1969 //
1970 //is in 10 slice?
1971 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1972 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1973 if (slice<0) slice=0;
1974 if (slice>10) slice=10;
1975 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1976 if (isOK) {
1977 fCurrentSlice=slice;
1978 fClustersCs = fClusters10[fCurrentSlice];
1979 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1980 fYcs = fY10[fCurrentSlice];
1981 fZcs = fZ10[fCurrentSlice];
1982 fNcs = fN10[fCurrentSlice];
1983 }
1984 }
1985 //
1986 //is in 5 slice?
1987 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1988 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1989 if (slice<0) slice=0;
1990 if (slice>5) slice=5;
1991 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
ae00569a 1992 if (isOK) {
e43c066c 1993 fCurrentSlice=slice;
1994 fClustersCs = fClusters5[fCurrentSlice];
1995 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1996 fYcs = fY5[fCurrentSlice];
1997 fZcs = fZ5[fCurrentSlice];
1998 fNcs = fN5[fCurrentSlice];
1999 }
2000 }
628e7bb0 2001 //
d9ead1a0 2002 fI = FindClusterIndex(fZmin);
2003 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2004 fSkip = 0;
2005 fAccepted = 0;
6518a6c5 2006
2007 return;
e43c066c 2008}
44347160 2009//------------------------------------------------------------------------
628e7bb0 2010Int_t AliITStrackerMI::AliITSlayer::
2011FindDetectorIndex(Double_t phi, Double_t z) const {
2012 //--------------------------------------------------------------------
2013 //This function finds the detector crossed by the track
2014 //--------------------------------------------------------------------
31fb8575 2015 Double_t dphi;
2016 if (fZOffset<0) // old geometry
2017 dphi = -(phi-fPhiOffset);
2018 else // new geometry
2019 dphi = phi-fPhiOffset;
2020
6518a6c5 2021
628e7bb0 2022 if (dphi < 0) dphi += 2*TMath::Pi();
2023 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2024 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2025 if (np>=fNladders) np-=fNladders;
2026 if (np<0) np+=fNladders;
2027
6518a6c5 2028
628e7bb0 2029 Double_t dz=fZOffset-z;
6518a6c5 2030 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2031 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
bf1cd124 2032 if (nz>=fNdetectors || nz<0) {
2033 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2034 return -1;
2035 }
628e7bb0 2036
23197852 2037 // ad hoc correction for 3rd ladder of SDD inner layer,
2038 // which is reversed (rotated by pi around local y)
2039 // this correction is OK only from AliITSv11Hybrid onwards
2040 if (GetR()>12. && GetR()<20.) { // SDD inner
2041 if(np==2) { // 3rd ladder
6801b453 2042 Double_t posMod252[3];
2043 AliITSgeomTGeo::GetTranslation(252,posMod252);
2044 // check the Z coordinate of Mod 252: if negative
2045 // (old SDD geometry in AliITSv11Hybrid)
2046 // the swap of numeration whould be applied
2047 if(posMod252[2]<0.){
2048 nz = (fNdetectors-1) - nz;
2049 }
2050 }
23197852 2051 }
2052 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2053
2054
628e7bb0 2055 return np*fNdetectors + nz;
2056}
44347160 2057//------------------------------------------------------------------------
ae00569a 2058const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2059{
e43c066c 2060 //--------------------------------------------------------------------
2061 // This function returns clusters within the "window"
2062 //--------------------------------------------------------------------
2063
ae00569a 2064 if (fCurrentSlice<0) {
e43c066c 2065 Double_t rpi2 = 2.*fR*TMath::Pi();
2066 for (Int_t i=fI; i<fImax; i++) {
2067 Double_t y = fY[i];
d9ead1a0 2068 Double_t z = fZ[i];
e43c066c 2069 if (fYmax<y) y -= rpi2;
1d4090b7 2070 if (fYmin>y) y += rpi2;
e43c066c 2071 if (y<fYmin) continue;
2072 if (y>fYmax) continue;
d9ead1a0 2073 // AD
2074 // skip clusters that are in "extended" road but they
2075 // 3sigma error does not touch the original road
2076 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2077 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2078 //
88e3c173 2079 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
e43c066c 2080 ci=i;
ae00569a 2081 if (!test) fI=i+1;
e43c066c 2082 return fClusters[i];
2083 }
ae00569a 2084 } else {
e43c066c 2085 for (Int_t i=fI; i<fImax; i++) {
2086 if (fYcs[i]<fYmin) continue;
2087 if (fYcs[i]>fYmax) continue;
88e3c173 2088 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
e43c066c 2089 ci=fClusterIndexCs[i];
ae00569a 2090 if (!test) fI=i+1;
e43c066c 2091 return fClustersCs[i];
2092 }
2093 }
2094 return 0;
2095}
44347160 2096//------------------------------------------------------------------------
e43c066c 2097Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2098const {
2099 //--------------------------------------------------------------------
afd25725 2100 // This function returns the layer thickness at this point (units X0)
e43c066c 2101 //--------------------------------------------------------------------
2102 Double_t d=0.0085;
e50912db 2103 x0=AliITSRecoParam::GetX0Air();
e43c066c 2104 if (43<fR&&fR<45) { //SSD2
2105 Double_t dd=0.0034;
2106 d=dd;
2107 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2108 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2109 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2110 for (Int_t i=0; i<12; i++) {
2111 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2112 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2113 d+=0.0034;
2114 break;
2115 }
2116 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2117 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2118 d+=0.0034;
2119 break;
2120 }
2121 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2122 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2123 }
2124 } else
2125 if (37<fR&&fR<41) { //SSD1
2126 Double_t dd=0.0034;
2127 d=dd;
2128 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2129 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2130 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2131 for (Int_t i=0; i<11; i++) {
2132 if (TMath::Abs(z-3.9*i)<0.15) {
2133 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2134 d+=dd;
2135 break;
2136 }
2137 if (TMath::Abs(z+3.9*i)<0.15) {
2138 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2139 d+=dd;
2140 break;
2141 }
2142 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2143 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2144 }
2145 } else
2146 if (13<fR&&fR<26) { //SDD
2147 Double_t dd=0.0033;
2148 d=dd;
2149 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2150
2151 if (TMath::Abs(y-1.80)<0.55) {
2152 d+=0.016;
2153 for (Int_t j=0; j<20; j++) {
2154 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2155 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2156 }
2157 }
2158 if (TMath::Abs(y+1.80)<0.55) {
2159 d+=0.016;
2160 for (Int_t j=0; j<20; j++) {
2161 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2162 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2163 }
2164 }
2165
2166 for (Int_t i=0; i<4; i++) {
2167 if (TMath::Abs(z-7.3*i)<0.60) {
2168 d+=dd;
2169 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2170 break;
2171 }
2172 if (TMath::Abs(z+7.3*i)<0.60) {
2173 d+=dd;
2174 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2175 break;
2176 }
2177 }
2178 } else
2179 if (6<fR&&fR<8) { //SPD2
2180 Double_t dd=0.0063; x0=21.5;
2181 d=dd;
2182 if (TMath::Abs(y-3.08)>0.5) d+=dd;
afd25725 2183 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
e43c066c 2184 } else
2185 if (3<fR&&fR<5) { //SPD1
2186 Double_t dd=0.0063; x0=21.5;
2187 d=dd;
2188 if (TMath::Abs(y+0.21)>0.6) d+=dd;
afd25725 2189 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
e43c066c 2190 }
2191
2192 return d;
2193}
44347160 2194//------------------------------------------------------------------------
23197852 2195AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2196fR(det.fR),
1c97ce2f 2197fRmisal(det.fRmisal),
23197852 2198fPhi(det.fPhi),
2199fSinPhi(det.fSinPhi),
2200fCosPhi(det.fCosPhi),
2201fYmin(det.fYmin),
2202fYmax(det.fYmax),
2203fZmin(det.fZmin),
2204fZmax(det.fZmax),
2205fIsBad(det.fIsBad),
2206fNChips(det.fNChips),
2207fChipIsBad(det.fChipIsBad)
2208{
2209 //Copy constructor
2210}
2211//------------------------------------------------------------------------
2212void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
a70ed6ad 2213 const AliITSDetTypeRec *detTypeRec)
23197852 2214{
2215 //--------------------------------------------------------------------
2216 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2217 //--------------------------------------------------------------------
2218
2219 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2220 // while in the tracker they start from 0 for each layer
2221 for(Int_t il=0; il<ilayer; il++)
2222 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2223
2224 Int_t detType;
2225 if (ilayer==0 || ilayer==1) { // ---------- SPD
2226 detType = 0;
2227 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2228 detType = 1;
2229 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2230 detType = 2;
2231 } else {
2232 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2233 return;
2234 }
2235
2236 // Get calibration from AliITSDetTypeRec
2237 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
d1181f10 2238 calib->SetModuleIndex(idet);
23197852 2239 AliITSCalibration *calibSPDdead = 0;
2240 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2241 if (calib->IsBad() ||
2242 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2243 {
2244 SetBad();
359ac74d 2245 // printf("lay %d bad %d\n",ilayer,idet);
23197852 2246 }
2247
2248 // Get segmentation from AliITSDetTypeRec
2249 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2250
2251 // Read info about bad chips
2252 fNChips = segm->GetMaximumChipIndex()+1;
2253 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2254 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2255 fChipIsBad = new Bool_t[fNChips];
2256 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2257 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2258 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
4187a601 2259 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
23197852 2260 }
2261
2262 return;
2263}
2264//------------------------------------------------------------------------
e50912db 2265Double_t AliITStrackerMI::GetEffectiveThickness()
e43c066c 2266{
2267 //--------------------------------------------------------------------
afd25725 2268 // Returns the thickness between the current layer and the vertex (units X0)
e43c066c 2269 //--------------------------------------------------------------------
e43c066c 2270
e50912db 2271 if(fUseTGeo!=0) {
2272 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2273 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2274 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2275 }
2276
afd25725 2277 // beam pipe
e50912db 2278 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2279 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
afd25725 2280
2281 // layers
2282 Double_t x0=0;
e43c066c 2283 Double_t xn=fgLayers[fI].GetR();
2284 for (Int_t i=0; i<fI; i++) {
2285 Double_t xi=fgLayers[i].GetR();
e50912db 2286 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2287 d+=dLayer*xi*xi;
e43c066c 2288 }
2289
afd25725 2290 // shields
e43c066c 2291 if (fI>1) {
e50912db 2292 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2293 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
e43c066c 2294 }
e43c066c 2295 if (fI>3) {
e50912db 2296 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2297 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
e43c066c 2298 }
e43c066c 2299 return d/(xn*xn);
2300}
44347160 2301//------------------------------------------------------------------------
e43c066c 2302Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
ae00569a 2303 //-------------------------------------------------------------------
e43c066c 2304 // This function returns number of clusters within the "window"
2305 //--------------------------------------------------------------------
2306 Int_t ncl=0;
2307 for (Int_t i=fI; i<fN; i++) {
00a7cc50 2308 const AliITSRecPoint *c=fClusters[i];
e43c066c 2309 if (c->GetZ() > fZmax) break;
2310 if (c->IsUsed()) continue;
2311 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2312 Double_t y=fR*det.GetPhi() + c->GetY();
2313
2314 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2315 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2316
2317 if (y<fYmin) continue;
2318 if (y>fYmax) continue;
2319 ncl++;
2320 }
2321 return ncl;
2322}
44347160 2323//------------------------------------------------------------------------
ae00569a 2324Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
7167ae53 2325 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
ae00569a 2326{
e43c066c 2327 //--------------------------------------------------------------------
ae00569a 2328 // This function refits the track "track" at the position "x" using
2329 // the clusters from "clusters"
ef7253ac 2330 // If "extra"==kTRUE,
ae00569a 2331 // the clusters from overlapped modules get attached to "track"
7167ae53 2332 // If "planeff"==kTRUE,
2333 // special approach for plane efficiency evaluation is applyed
e43c066c 2334 //--------------------------------------------------------------------
6518a6c5 2335
e50912db 2336 Int_t index[AliITSgeomTGeo::kNLayers];
e43c066c 2337 Int_t k;
e50912db 2338 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
ae00569a 2339 Int_t nc=clusters->GetNumberOfClusters();
e43c066c 2340 for (k=0; k<nc; k++) {
ae00569a 2341 Int_t idx=clusters->GetClusterIndex(k);
2342 Int_t ilayer=(idx&0xf0000000)>>28;
2343 index[ilayer]=idx;
e43c066c 2344 }
67c1e979 2345
7167ae53 2346 return RefitAt(xx,track,index,extra,planeeff); // call the method below
67c1e979 2347}
44347160 2348//------------------------------------------------------------------------
ae00569a 2349Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
7167ae53 2350 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
ae00569a 2351{
67c1e979 2352 //--------------------------------------------------------------------
ae00569a 2353 // This function refits the track "track" at the position "x" using
67c1e979 2354 // the clusters from array
ae00569a 2355 // If "extra"==kTRUE,
2356 // the clusters from overlapped modules get attached to "track"
7167ae53 2357 // If "planeff"==kTRUE,
2358 // special approach for plane efficiency evaluation is applyed
67c1e979 2359 //--------------------------------------------------------------------
e50912db 2360 Int_t index[AliITSgeomTGeo::kNLayers];
67c1e979 2361 Int_t k;
e50912db 2362 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
67c1e979 2363 //
e50912db 2364 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
ae00569a 2365 index[k]=clusters[k];
67c1e979 2366 }
2367
74008278 2368 // special for cosmics and TPC prolonged tracks:
269709f8 2369 // propagate to the innermost of:
2370 // - innermost layer crossed by the track
2371 // - innermost layer where a cluster was associated to the track
74008278 2372 static AliITSRecoParam *repa = NULL;
2373 if(!repa){
2374 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2375 if(!repa){
2376 repa = AliITSRecoParam::GetHighFluxParam();
2377 AliWarning("Using default AliITSRecoParam class");
2378 }
2379 }
2380 Int_t evsp=repa->GetEventSpecie();
15f91fcd 2381 ULong_t trStatus=0;
2382 if(track->GetESDtrack()) trStatus=track->GetStatus();
74008278 2383 Int_t innermostlayer=0;
15f91fcd 2384 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
74008278 2385 innermostlayer=5;
2386 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2387 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
269709f8 2388 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2389 index[innermostlayer] >= 0 ) break;
74008278 2390 }
269709f8 2391
74008278 2392 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
ae00569a 2393 }
afd25725 2394
ae00569a 2395 Int_t modstatus=1; // found
2396 Float_t xloc,zloc;
67c1e979 2397 Int_t from, to, step;
ae00569a 2398 if (xx > track->GetX()) {
e50912db 2399 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
67c1e979 2400 step=+1;
2401 } else {
e50912db 2402 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
67c1e979 2403 step=-1;
2404 }
ae00569a 2405 TString dir = (step>0 ? "outward" : "inward");
67c1e979 2406
ae00569a 2407 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
6518a6c5 2408 AliITSlayer &layer=fgLayers[ilayer];
67c1e979 2409 Double_t r=layer.GetR();
4fa7f7d1 2410
afd25725 2411 if (step<0 && xx>r) break;
2412
2413 // material between SSD and SDD, SDD and SPD
6518a6c5 2414 Double_t hI=ilayer-0.5*step;
e50912db 2415 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
ae00569a 2416 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
e50912db 2417 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
ae00569a 2418 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
67c1e979 2419
ae00569a 2420
afd25725 2421 Double_t oldGlobXYZ[3];
d1181f10 2422 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
afd25725 2423
67178f34 2424 // continue if we are already beyond this layer
2425 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2426 if(step>0 && oldGlobR > r) continue; // going outward
2427 if(step<0 && oldGlobR < r) continue; // going inward
2428
8602c008 2429 Double_t phi,z;
ae00569a 2430 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
67c1e979 2431
67c1e979 2432 Int_t idet=layer.FindDetectorIndex(phi,z);
6518a6c5 2433
2434 // check if we allow a prolongation without point for large-eta tracks
ae00569a 2435 Int_t skip = CheckSkipLayer(track,ilayer,idet);
6518a6c5 2436 if (skip==2) {
ae00569a 2437 modstatus = 4; // out in z
d1181f10 2438 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2439 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2440 }
7131593b 2441 // cross layer
2442 // apply correction for material of the current layer
2443 // add time if going outward
2444 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
6518a6c5 2445 continue;
2446 }
2447
e50912db 2448 if (idet<0) return kFALSE;
ae00569a 2449
4fa7f7d1 2450
67c1e979 2451 const AliITSdetector &det=layer.GetDetector(idet);
545af9ef 2452 // only for ITS-SA tracks refit
2453 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2454 //
1c97ce2f 2455 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
23197852 2456
ae00569a 2457 track->SetDetectorIndex(idet);
d1181f10 2458 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
67c1e979 2459
23197852 2460 Double_t dz,zmin,zmax,dy,ymin,ymax;
ae00569a 2461
2462 const AliITSRecPoint *clAcc=0;
44347160 2463 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
67c1e979 2464
6518a6c5 2465 Int_t idx=index[ilayer];
ae00569a 2466 if (idx>=0) { // cluster in this layer
2467 modstatus = 6; // no refit
2468 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2469 if (cl) {
2470 if (idet != cl->GetDetectorIndex()) {
2471 idet=cl->GetDetectorIndex();
a5d0fae6 2472 const AliITSdetector &detc=layer.GetDetector(idet);
2473 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
ae00569a 2474 track->SetDetectorIndex(idet);
d1181f10 2475 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
ae00569a 2476 }
ae00569a 2477 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2478 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2479 if (chi2<maxchi2) {
2480 clAcc=cl;
2481 maxchi2=chi2;
2482 modstatus = 1; // found
2483 } else {
7167ae53 2484 return kFALSE; //
ae00569a 2485 }
2486 }
2487 } else { // no cluster in this layer
2488 if (skip==1) {
2489 modstatus = 3; // skipped
0ed58a47 2490 // Plane Eff determination:
2491 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2492 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2493 UseTrackForPlaneEff(track,ilayer);
7167ae53 2494 }
ae00569a 2495 } else {
2496 modstatus = 5; // no cls in road
2497 // check dead
1c97ce2f 2498 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2499 dz = 0.5*(zmax-zmin);
2500 dy = 0.5*(ymax-ymin);
2501 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
ae00569a 2502 if (dead==1) modstatus = 7; // holes in z in SPD
f8720bda 2503 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
ae00569a 2504 }
e43c066c 2505 }
ae00569a 2506
2507 if (clAcc) {
2508 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
09cf9d66 2509 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
e43c066c 2510 }
ae00569a 2511 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2512
2513
ef3508c5 2514 if (extra && clAcc) { // search for extra clusters in overlapped modules
ae00569a 2515 AliITStrackV2 tmp(*track);
1c97ce2f 2516 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2517 layer.SelectClusters(zmin,zmax,ymin,ymax);
ae00569a 2518
2519 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2520 Int_t idetExtra=-1;
a5d0fae6 2521 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2522 Double_t tolerance=0.1;
ae00569a 2523 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2524 // only clusters in another module! (overlaps)
2525 idetExtra = clExtra->GetDetectorIndex();
2526 if (idet == idetExtra) continue;
2527
a5d0fae6 2528 const AliITSdetector &detx=layer.GetDetector(idetExtra);
ae00569a 2529
1c97ce2f 2530 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
ae00569a 2531 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2532 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
1c97ce2f 2533 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
ae00569a 2534
2535 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2536 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2537 }
2538 if (cci>=0) {
2539 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2540 track->SetExtraModule(ilayer,idetExtra);
2541 }
2542 } // end search for extra clusters in overlapped modules
1c97ce2f 2543
afd25725 2544 // Correct for material of the current layer
7131593b 2545 // cross material
2546 // add time if going outward
ae00569a 2547 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
545af9ef 2548 track->SetCheckInvariant(kTRUE);
ae00569a 2549 } // end loop on layers
2550
db355ee7 2551 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
e43c066c 2552
e43c066c 2553 return kTRUE;
2554}
44347160 2555//------------------------------------------------------------------------
15dd636f 2556Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
e43c066c 2557{
2558 //
2559 // calculate normalized chi2
2560 // return NormalizedChi2(track,0);
2561 Float_t chi2 = 0;
2562 Float_t sum=0;
2563 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2564 // track->fdEdxMismatch=0;
2565 Float_t dedxmismatch =0;
2566 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2567 if (mode<100){
2568 for (Int_t i = 0;i<6;i++){
9fd412a8 2569 if (track->GetClIndex(i)>=0){
e43c066c 2570 Float_t cerry, cerrz;
2571 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2572 else
b9671574 2573 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
e43c066c 2574 cerry*=cerry;
2575 cerrz*=cerrz;
b9671574 2576 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
572f41f9 2577 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
b9671574 2578 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
e43c066c 2579 if (ratio<0.5) {
2580 cchi2+=(0.5-ratio)*10.;
2581 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2582 dedxmismatch+=(0.5-ratio)*10.;
2583 }
2584 }
2585 if (i<2 ||i>3){
b9671574 2586 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
e43c066c 2587 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2588 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2589 if (i<2) chi2+=2*cl->GetDeltaProbability();
2590 }
2591 chi2+=cchi2;
2592 sum++;
2593 }
2594 }
b9671574 2595 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2596 track->SetdEdxMismatch(dedxmismatch);
e43c066c 2597 }
2598 }
2599 else{
2600 for (Int_t i = 0;i<4;i++){
9fd412a8 2601 if (track->GetClIndex(i)>=0){
e43c066c 2602 Float_t cerry, cerrz;
2603 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
b9671574 2604 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
e43c066c 2605 cerry*=cerry;
2606 cerrz*=cerrz;
b9671574 2607 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2608 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
e43c066c 2609 sum++;
2610 }
2611 }
2612 for (Int_t i = 4;i<6;i++){
9fd412a8 2613 if (track->GetClIndex(i)>=0){
e43c066c 2614 Float_t cerry, cerrz;
2615 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
b9671574 2616 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
e43c066c 2617 cerry*=cerry;
2618 cerrz*=cerrz;
2619 Float_t cerryb, cerrzb;
2620 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
b9671574 2621 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
e43c066c 2622 cerryb*=cerryb;
2623 cerrzb*=cerrzb;
b9671574 2624 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2625 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
e43c066c 2626 sum++;
2627 }
2628 }
2629 }
b9671574 2630 if (track->GetESDtrack()->GetTPCsignal()>85){
2631 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
e43c066c 2632 if (ratio<0.5) {
2633 chi2+=(0.5-ratio)*5.;
2634 }
2635 if (ratio>2){
2636 chi2+=(ratio-2.0)*3;
2637 }
2638 }
2639 //
b9671574 2640 Double_t match = TMath::Sqrt(track->GetChi22());
2641 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2755f080 2642 if (!track->GetConstrain()) {
2643 if (track->GetNumberOfClusters()>2) {
2644 match/=track->GetNumberOfClusters()-2.;
2645 } else {
2646 match=0;
2647 }
2648 }
e43c066c 2649 if (match<0) match=0;
f8720bda 2650
2651 // penalty factor for missing points (NDeadZone>0), but no penalty
2652 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2653 // or there is a dead from OCDB)
2654 Float_t deadzonefactor = 0.;
2655 if (track->GetNDeadZone()>0.) {
bf1cd124 2656 Int_t sumDeadZoneProbability=0;
2657 for(Int_t ilay=0;ilay<6;ilay++) {
2658 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2659 }
2660 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2661 if(nDeadZoneWithProbNot1>0) {
2662 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
a5a317a9 2663 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
bf1cd124 2664 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2665 Float_t one = 1.;
2666 deadZoneProbability = TMath::Min(deadZoneProbability,one);
f8720bda 2667 deadzonefactor = 3.*(1.1-deadZoneProbability);
2668 }
2669 }
2670
b9671574 2671 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2672 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2673 1./(1.+track->GetNSkipped()));
a5a317a9 2674 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
bf1cd124 2675 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2676 return normchi2;
e43c066c 2677}
44347160 2678//------------------------------------------------------------------------
a70ed6ad 2679Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
e43c066c 2680{
2681 //
2682 // return matching chi2 between two tracks
d1181f10 2683 Double_t largeChi2=1000.;
2684
15dd636f 2685 AliITStrackMI track3(*track2);
d1181f10 2686 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
e43c066c 2687 TMatrixD vec(5,1);
6c94f330 2688 vec(0,0)=track1->GetY() - track3.GetY();
2689 vec(1,0)=track1->GetZ() - track3.GetZ();
2690 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2691 vec(3,0)=track1->GetTgl() - track3.GetTgl();
6c23ffed 2692 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
e43c066c 2693 //
2694 TMatrixD cov(5,5);
6c94f330 2695 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2696 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2697 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2698 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2699 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
e43c066c 2700
6c94f330 2701 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2702 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2703 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2704 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
e43c066c 2705 //
6c94f330 2706 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2707 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2708 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
e43c066c 2709 //
6c94f330 2710 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2711 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
e43c066c 2712 //
6c94f330 2713 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
e43c066c 2714
2715 cov.Invert();
2716 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2717 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2718 return chi2(0,0);
2719}
44347160 2720//------------------------------------------------------------------------
a70ed6ad 2721Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
e43c066c 2722{
2723 //
afd25725 2724 // return probability that given point (characterized by z position and error)
2725 // is in SPD dead zone
bf1cd124 2726 // This method assumes that fSPDdetzcentre is ordered from -z to +z
e43c066c 2727 //
afd25725 2728 Double_t probability = 0.;
bf1cd124 2729 Double_t nearestz = 0.,distz=0.;
2730 Int_t nearestzone = -1;
2731 Double_t mindistz = 1000.;
2732
2733 // find closest dead zone
2734 for (Int_t i=0; i<3; i++) {
2735 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2736 if (distz<mindistz) {
2737 nearestzone=i;
2738 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2739 mindistz=distz;
2740 }
2741 }
2742
2743 // too far from dead zone
2744 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2745
2746
afd25725 2747 Double_t zmin, zmax;
bf1cd124 2748 if (nearestzone==0) { // dead zone at z = -7
e50912db 2749 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2750 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
bf1cd124 2751 } else if (nearestzone==1) { // dead zone at z = 0
e50912db 2752 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2753 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
bf1cd124 2754 } else if (nearestzone==2) { // dead zone at z = +7
2755 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2756 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
afd25725 2757 } else {
2758 zmin = 0.;
2759 zmax = 0.;
e43c066c 2760 }
afd25725 2761 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2762 // dead zone)
bb7e41dd 2763 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2764 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
bf1cd124 2765 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
e43c066c 2766 return probability;
2767}
44347160 2768//------------------------------------------------------------------------
a70ed6ad 2769Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
e43c066c 2770{
2771 //
2772 // calculate normalized chi2
2773 Float_t chi2[6];
2774 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2775 Float_t ncl = 0;
2776 for (Int_t i = 0;i<6;i++){
b9671574 2777 if (TMath::Abs(track->GetDy(i))>0){
2778 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2779 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
e43c066c 2780 ncl++;
2781 }
2782 else{chi2[i]=10000;}
2783 }
2784 Int_t index[6];
2785 TMath::Sort(6,chi2,index,kFALSE);
2786 Float_t max = float(ncl)*fac-1.;
2787 Float_t sumchi=0, sumweight=0;
2788 for (Int_t i=0;i<max+1;i++){
2789 Float_t weight = (i<max)?1.:(max+1.-i);
2790 sumchi+=weight*chi2[index[i]];
2791 sumweight+=weight;
2792 }
2793 Double_t normchi2 = sumchi/sumweight;
2794 return normchi2;
2795}
44347160 2796//------------------------------------------------------------------------
a70ed6ad 2797Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
e43c066c 2798{
2799 //
2800 // calculate normalized chi2
2801 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2802 Int_t npoints = 0;
2803 Double_t res =0;
2804 for (Int_t i=0;i<6;i++){
b9671574 2805 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2806 Double_t sy1 = forwardtrack->GetSigmaY(i);
2807 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2808 Double_t sy2 = backtrack->GetSigmaY(i);
2809 Double_t sz2 = backtrack->GetSigmaZ(i);
e43c066c 2810 if (i<2){ sy2=1000.;sz2=1000;}
2811 //
b9671574 2812 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2813 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
e43c066c 2814 //
2815 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2816 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2817 //
2818 res+= nz0*nz0+ny0*ny0;
2819 npoints++;
2820 }
2821 if (npoints>1) return
6c23ffed 2822 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
e43c066c 2823 //2*forwardtrack->fNUsed+
b9671574 2824 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2825 1./(1.+forwardtrack->GetNSkipped()));
e43c066c 2826 return 1000;
2827}
44347160 2828//------------------------------------------------------------------------
e43c066c 2829Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2830 //--------------------------------------------------------------------
2831 // Return pointer to a given cluster
2832 //--------------------------------------------------------------------
2833 Int_t l=(index & 0xf0000000) >> 28;
2834 Int_t c=(index & 0x0fffffff) >> 00;
2835 return fgLayers[l].GetWeight(c);
2836}
44347160 2837//------------------------------------------------------------------------
a70ed6ad 2838void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
e43c066c 2839{
2840 //---------------------------------------------
2841 // register track to the list
628e7bb0 2842 //
b9671574 2843 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
628e7bb0 2844 //
2845 //
e43c066c 2846 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2847 Int_t index = track->GetClusterIndex(icluster);
2848 Int_t l=(index & 0xf0000000) >> 28;
2849 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2850 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2851 for (Int_t itrack=0;itrack<4;itrack++){
b9671574 2852 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2853 fgLayers[l].SetClusterTracks(itrack,c,id);
e43c066c 2854 break;
2855 }
2856 }
2857 }
2858}
44347160 2859//------------------------------------------------------------------------
a70ed6ad 2860void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
e43c066c 2861{
2862 //---------------------------------------------
2863 // unregister track from the list
2864 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2865 Int_t index = track->GetClusterIndex(icluster);
2866 Int_t l=(index & 0xf0000000) >> 28;
2867 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2868 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2869 for (Int_t itrack=0;itrack<4;itrack++){
b9671574 2870 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2871 fgLayers[l].SetClusterTracks(itrack,c,-1);
e43c066c 2872 }
2873 }
2874 }
2875}
44347160 2876//------------------------------------------------------------------------
00a7cc50 2877Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
e43c066c 2878{
2879 //-------------------------------------------------------------
2880 //get number of shared clusters
2881 //-------------------------------------------------------------
2882 Float_t shared=0;
2883 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2884 // mean number of clusters
2885 Float_t *ny = GetNy(id), *nz = GetNz(id);
2886
2887
2888 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2889 Int_t index = track->GetClusterIndex(icluster);
2890 Int_t l=(index & 0xf0000000) >> 28;
2891 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 2892 if (c>fgLayers[l].GetNumberOfClusters()) continue;
5668b5b2 2893 // if (ny[l]<1.e-13){
2894 // printf("problem\n");
2895 // }
00a7cc50 2896 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
e43c066c 2897 Float_t weight=1;
2898 //
2899 Float_t deltan = 0;
2900 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
572f41f9 2901 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2902 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
e43c066c 2903 if (l<2 || l>3){
2904 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2905 }
2906 else{
2907 deltan = (cl->GetNz()-nz[l]);
2908 }
2909 if (deltan>2.0) continue; // extended - highly probable shared cluster
2910 weight = 2./TMath::Max(3.+deltan,2.);
2911 //
2912 for (Int_t itrack=0;itrack<4;itrack++){
b9671574 2913 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
e43c066c 2914 list[l]=index;
00a7cc50 2915 clist[l] = (AliITSRecPoint*)GetCluster(index);
25015f7a 2916 track->SetSharedWeight(l,weight);
e43c066c 2917 shared+=weight;
2918 break;
2919 }
2920 }
2921 }
b9671574 2922 track->SetNUsed(shared);
e43c066c 2923 return shared;
2924}
44347160 2925//------------------------------------------------------------------------
a70ed6ad 2926Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
e43c066c 2927{
2928 //
2929 // find first shared track
2930 //
2931 // mean number of clusters
2932 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2933 //
2934 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2935 Int_t sharedtrack=100000;
2936 Int_t tracks[24],trackindex=0;
2937 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2938 //
2939 for (Int_t icluster=0;icluster<6;icluster++){
2940 if (clusterlist[icluster]<0) continue;
2941 Int_t index = clusterlist[icluster];
2942 Int_t l=(index & 0xf0000000) >> 28;
2943 Int_t c=(index & 0x0fffffff) >> 00;
5668b5b2 2944 // if (ny[l]<1.e-13){
2945 // printf("problem\n");
2946 // }
b9671574 2947 if (c>fgLayers[l].GetNumberOfClusters()) continue;
e43c066c 2948 //if (l>3) continue;
00a7cc50 2949 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
e43c066c 2950 //
2951 Float_t deltan = 0;
2952 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
572f41f9 2953 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2954 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
e43c066c 2955 if (l<2 || l>3){
2956 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2957 }
2958 else{
2959 deltan = (cl->GetNz()-nz[l]);
2960 }
2961 if (deltan>2.0) continue; // extended - highly probable shared cluster
2962 //
2963 for (Int_t itrack=3;itrack>=0;itrack--){
b9671574 2964 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2965 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2966 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
e43c066c 2967 trackindex++;
2968 }
2969 }
2970 }
2971 if (trackindex==0) return -1;
2972 if (trackindex==1){
2973 sharedtrack = tracks[0];
2974 }else{
2975 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2976 else{
2977 //
a5d0fae6 2978 Int_t tracks2[24], cluster[24];
2979 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
e43c066c 2980 Int_t index =0;
2981 //
2982 for (Int_t i=0;i<trackindex;i++){
2983 if (tracks[i]<0) continue;
a5d0fae6 2984 tracks2[index] = tracks[i];
e43c066c 2985 cluster[index]++;
2986 for (Int_t j=i+1;j<trackindex;j++){
2987 if (tracks[j]<0) continue;
2988 if (tracks[j]==tracks[i]){
2989 cluster[index]++;
2990 tracks[j]=-1;
2991 }
2992 }
2993 index++;
2994 }
2995 Int_t max=0;
2996 for (Int_t i=0;i<index;i++){
2997 if (cluster[index]>max) {
a5d0fae6 2998 sharedtrack=tracks2[index];
e43c066c 2999 max=cluster[index];
3000 }
3001 }
3002 }
3003 }
3004
3005 if (sharedtrack>=100000) return -1;
3006 //
3007 // make list of overlaps
3008 shared =0;
3009 for (Int_t icluster=0;icluster<6;icluster++){
3010 if (clusterlist[icluster]<0) continue;
3011 Int_t index = clusterlist[icluster];
3012 Int_t l=(index & 0xf0000000) >> 28;
3013 Int_t c=(index & 0x0fffffff) >> 00;
b9671574 3014 if (c>fgLayers[l].GetNumberOfClusters()) continue;
00a7cc50 3015 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
e43c066c 3016 if (l==0 || l==1){
3017 if (cl->GetNy()>2) continue;
3018 if (cl->GetNz()>2) continue;
3019 }
3020 if (l==4 || l==5){
3021 if (cl->GetNy()>3) continue;
3022 if (cl->GetNz()>3) continue;
3023 }
3024 //
3025 for (Int_t itrack=3;itrack>=0;itrack--){
b9671574 3026 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3027 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
e43c066c 3028 overlist[l]=index;
3029 shared++;
3030 }
3031 }
3032 }
3033 return sharedtrack;
3034}
44347160 3035//------------------------------------------------------------------------
15dd636f 3036AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
e43c066c 3037 //
3038 // try to find track hypothesys without conflicts
3039 // with minimal chi2;
3040 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3041 Int_t entries1 = arr1->GetEntriesFast();
3042 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
15dd636f 3043 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
e43c066c 3044 Int_t entries2 = arr2->GetEntriesFast();
15dd636f 3045 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
e43c066c 3046 //
15dd636f 3047 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3048 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
6c23ffed 3049 if (track10->Pt()>0.5+track20->Pt()) return track10;
e43c066c 3050
3051 for (Int_t itrack=0;itrack<entries1;itrack++){
15dd636f 3052 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
e43c066c 3053 UnRegisterClusterTracks(track,trackID1);
3054 }
3055 //
3056 for (Int_t itrack=0;itrack<entries2;itrack++){
15dd636f 3057 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
e43c066c 3058 UnRegisterClusterTracks(track,trackID2);
3059 }
3060 Int_t index1=0;
3061 Int_t index2=0;
3062 Float_t maxconflicts=6;
3063 Double_t maxchi2 =1000.;
3064 //
3065 // get weight of hypothesys - using best hypothesy
3066 Double_t w1,w2;
3067
3068 Int_t list1[6],list2[6];
00a7cc50 3069 AliITSRecPoint *clist1[6], *clist2[6] ;
e43c066c 3070 RegisterClusterTracks(track10,trackID1);
3071 RegisterClusterTracks(track20,trackID2);
3072 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3073 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3074 UnRegisterClusterTracks(track10,trackID1);
3075 UnRegisterClusterTracks(track20,trackID2);
3076 //
3077 // normalized chi2
3078 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3079 Float_t nerry[6],nerrz[6];
3080 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3081 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3082 for (Int_t i=0;i<6;i++){
3083 if ( (erry1[i]>0) && (erry2[i]>0)) {
3084 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3085 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3086 }else{
3087 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3088 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3089 }
b9671574 3090 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3091 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3092 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
e43c066c 3093 ncl1++;
3094 }
b9671574 3095 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3096 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3097 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
e43c066c 3098 ncl2++;
3099 }
3100 }
3101 chi21/=ncl1;
3102 chi22/=ncl2;
3103 //
3104 //
b9671574 3105 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3106 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
e43c066c 3107 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3108 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3109 //
3110 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3111 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
6c23ffed 3112 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
e43c066c 3113 );
3114 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3115 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
6c23ffed 3116 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
e43c066c 3117 );
3118
3119 Double_t sumw = w1+w2;
3120 w1/=sumw;
3121 w2/=sumw;
3122 if (w1<w2*0.5) {
3123 w1 = (d2+0.5)/(d1+d2+1.);
3124 w2 = (d1+0.5)/(d1+d2+1.);
3125 }
3126 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3127 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3128 //
3129 // get pair of "best" hypothesys
3130 //
3131 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3132 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3133
3134 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
15dd636f 3135 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
e43c066c 3136 //if (track1->fFakeRatio>0) continue;
3137 RegisterClusterTracks(track1,trackID1);
3138 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
15dd636f 3139 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
e43c066c 3140
3141 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3142 //if (track2->fFakeRatio>0) continue;
3143 Float_t nskipped=0;
3144 RegisterClusterTracks(track2,trackID2);
e43c066c 3145 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3146 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3147 UnRegisterClusterTracks(track2,trackID2);
3148 //
b9671574 3149 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3150 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
e43c066c 3151 if (nskipped>0.5) continue;
3152 //
3153 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3154 if (conflict1+1<cconflict1) continue;
3155 if (conflict2+1<cconflict2) continue;
3156 Float_t conflict=0;
3157 Float_t sumchi2=0;
3158 Float_t sum=0;
3159 for (Int_t i=0;i<6;i++){
3160 //
3161 Float_t c1 =0.; // conflict coeficients
3162 Float_t c2 =0.;
3163 if (clist1[i]&&clist2[i]){
3164 Float_t deltan = 0;
3165 if (i<2 || i>3){
3166 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3167 }
3168 else{
3169 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3170 }
3171 c1 = 2./TMath::Max(3.+deltan,2.);
3172 c2 = 2./TMath::Max(3.+deltan,2.);
3173 }
3174 else{
3175 if (clist1[i]){
3176 Float_t deltan = 0;
3177 if (i<2 || i>3){
3178 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3179 }
3180 else{
3181 deltan = (clist1[i]->GetNz()-nz1[i]);
3182 }
3183 c1 = 2./TMath::Max(3.+deltan,2.);
3184 c2 = 0;
3185 }
3186
3187 if (clist2[i]){
3188 Float_t deltan = 0;
3189 if (i<2 || i>3){
3190 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3191 }
3192 else{
3193 deltan = (clist2[i]->GetNz()-nz2[i]);
3194 }
3195 c2 = 2./TMath::Max(3.+deltan,2.);
3196 c1 = 0;
3197 }
3198 }
3199 //
a5d0fae6 3200 chi21=0;chi22=0;
b9671574 3201 if (TMath::Abs(track1->GetDy(i))>0.) {
3202 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3203 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
e43c066c 3204 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
b9671574 3205 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
e43c066c 3206 }else{
b9671574 3207 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
e43c066c 3208 }
3209 //
b9671574 3210 if (TMath::Abs(track2->GetDy(i))>0.) {
3211 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3212 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
e43c066c 3213 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3214 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3215 }
3216 else{
b9671574 3217 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
e43c066c 3218 }
3219 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3220 if (chi21>0) sum+=w1;
3221 if (chi22>0) sum+=w2;
3222 conflict+=(c1+c2);
3223 }
b9671574 3224 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3225 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
e43c066c 3226 Double_t normchi2 = 2*conflict+sumchi2/sum;
3227 if ( normchi2 <maxchi2 ){
3228 index1 = itrack1;
3229 index2 = itrack2;
3230 maxconflicts = conflict;
3231 maxchi2 = normchi2;
3232 }
3233 }
3234 UnRegisterClusterTracks(track1,trackID1);
3235 }
3236 //
3237 // if (maxconflicts<4 && maxchi2<th0){
3238 if (maxchi2<th0*2.){
b9671574 3239 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
15dd636f 3240 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
b9671574 3241 track1->SetChi2MIP(5,maxconflicts);
3242 track1->SetChi2MIP(6,maxchi2);
3243 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
e43c066c 3244 // track1->UpdateESDtrack(AliESDtrack::kITSin);
b9671574 3245 track1->SetChi2MIP(8,index1);
e43c066c 3246 fBestTrackIndex[trackID1] =index1;
3247 UpdateESDtrack(track1, AliESDtrack::kITSin);
3248 }
b9671574 3249 else if (track10->GetChi2MIP(0)<th1){
3250 track10->SetChi2MIP(5,maxconflicts);
3251 track10->SetChi2MIP(6,maxchi2);
e43c066c 3252 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3253 UpdateESDtrack(track10,AliESDtrack::kITSin);
3254 }
3255
3256 for (Int_t itrack=0;itrack<entries1;itrack++){
15dd636f 3257 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
e43c066c 3258 UnRegisterClusterTracks(track,trackID1);
3259 }
3260 //
3261 for (Int_t itrack=0;itrack<entries2;itrack++){
15dd636f 3262 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
e43c066c 3263 UnRegisterClusterTracks(track,trackID2);
3264 }
3265
44347160 3266 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3267 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3268 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3269 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
e43c066c 3270 RegisterClusterTracks(track10,trackID1);
3271 }
44347160 3272 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3273 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3274 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3275 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
e43c066c 3276 RegisterClusterTracks(track20,trackID2);
3277 }
3278 return track10;
3279
3280}
44347160 3281//------------------------------------------------------------------------
e43c066c 3282void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3283 //--------------------------------------------------------------------
3284 // This function marks clusters assigned to the track
3285 //--------------------------------------------------------------------
3286 AliTracker::UseClusters(t,from);
3287
00a7cc50 3288 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
e43c066c 3289 //if (c->GetQ()>2) c->Use();
3290 if (c->GetSigmaZ2()>0.1) c->Use();
00a7cc50 3291 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
e43c066c 3292 //if (c->GetQ()>2) c->Use();
3293 if (c->GetSigmaZ2()>0.1) c->Use();
3294
3295}
44347160 3296//------------------------------------------------------------------------
15dd636f 3297void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
e43c066c 3298{
3299 //------------------------------------------------------------------
3300 // add track to the list of hypothesys
3301 //------------------------------------------------------------------
3302
fddf8459 3303 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3304 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
e43c066c 3305 //
3306 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3307 if (!array) {
3308 array = new TObjArray(10);
3309 fTrackHypothesys.AddAt(array,esdindex);
3310 }
3311 array->AddLast(track);
3312}
44347160 3313//------------------------------------------------------------------------
e43c066c 3314void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3315{
3316 //-------------------------------------------------------------------
3317 // compress array of track hypothesys
3318 // keep only maxsize best hypothesys
3319 //-------------------------------------------------------------------
3320 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3321 if (! (fTrackHypothesys.At(esdindex)) ) return;
3322 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3323 Int_t entries = array->GetEntriesFast();
3324 //
3325 //- find preliminary besttrack as a reference
3326 Float_t minchi2=10000;
3327 Int_t maxn=0;
15dd636f 3328 AliITStrackMI * besttrack=0;
e43c066c 3329 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
15dd636f 3330 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
e43c066c 3331 if (!track) continue;
3332 Float_t chi2 = NormalizedChi2(track,0);
3333 //
b9671574 3334 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
e43c066c 3335 track->SetLabel(tpcLabel);
3336 CookdEdx(track);
b9671574 3337 track->SetFakeRatio(1.);
e43c066c 3338 CookLabel(track,0.); //For comparison only
3339 //
44347160 3340 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3341 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
e43c066c 3342 if (track->GetNumberOfClusters()<maxn) continue;
3343 maxn = track->GetNumberOfClusters();
3344 if (chi2<minchi2){
3345 minchi2=chi2;
3346 besttrack=track;
3347 }
3348 }
3349 else{
b9671574 3350 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 3351 delete array->RemoveAt(itrack);
3352 }
e43c066c 3353 }
3354 }
3355 if (!besttrack) return;
3356 //
3357 //
3358 //take errors of best track as a reference
3359 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3360 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
a5d0fae6 3361 for (Int_t j=0;j<6;j++) {
9fd412a8 3362 if (besttrack->GetClIndex(j)>=0){
a5d0fae6 3363 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3364 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3365 ny[j] = besttrack->GetNy(j);
3366 nz[j] = besttrack->GetNz(j);
e43c066c 3367 }
3368 }
3369 //
3370 // calculate normalized chi2
3371 //
3372 Float_t * chi2 = new Float_t[entries];
3373 Int_t * index = new Int_t[entries];
3374 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3375 for (Int_t itrack=0;itrack<entries;itrack++){
15dd636f 3376 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
e43c066c 3377 if (track){
bf1cd124 3378 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
b9671574 3379 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
44347160 3380 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
b9671574 3381 chi2[itrack] = track->GetChi2MIP(0);
628e7bb0 3382 else{
b9671574 3383 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 3384 delete array->RemoveAt(itrack);
3385 }
3386 }
e43c066c 3387 }
3388 }
3389 //
3390 TMath::Sort(entries,chi2,index,kFALSE);
15dd636f 3391 besttrack = (AliITStrackMI*)array->At(index[0]);
bf1cd124 3392 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
44347160 3393 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
a5d0fae6 3394 for (Int_t j=0;j<6;j++){
9fd412a8 3395 if (besttrack->GetClIndex(j)>=0){
a5d0fae6 3396 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3397 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3398 ny[j] = besttrack->GetNy(j);
3399 nz[j] = besttrack->GetNz(j);
e43c066c 3400 }
3401 }
3402 }
3403 //
3404 // calculate one more time with updated normalized errors
3405 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3406 for (Int_t itrack=0;itrack<entries;itrack++){
15dd636f 3407 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
e43c066c 3408 if (track){
bf1cd124 3409 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3410 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
44347160 3411 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
b9671574 3412 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
628e7bb0 3413 else
3414 {
b9671574 3415 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 3416 delete array->RemoveAt(itrack);
3417 }
3418 }
e43c066c 3419 }
3420 }
3421 entries = array->GetEntriesFast();
3422 //
628e7bb0 3423 //
e43c066c 3424 if (entries>0){
3425 TObjArray * newarray = new TObjArray();
3426 TMath::Sort(entries,chi2,index,kFALSE);
15dd636f 3427 besttrack = (AliITStrackMI*)array->At(index[0]);
e43c066c 3428 if (besttrack){
bf1cd124 3429 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
e43c066c 3430 //
a5d0fae6 3431 for (Int_t j=0;j<6;j++){
3432 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3433 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3434 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3435 ny[j] = besttrack->GetNy(j);
3436 nz[j] = besttrack->GetNz(j);
e43c066c 3437 }
3438 }
b9671574 3439 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
a5d0fae6 3440 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
e43c066c 3441 Float_t minn = besttrack->GetNumberOfClusters()-3;
3442 Int_t accepted=0;
3443 for (Int_t i=0;i<entries;i++){
15dd636f 3444 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
e43c066c 3445 if (!track) continue;
3446 if (accepted>maxcut) break;
b9671574 3447 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3448 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3449 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
628e7bb0 3450 delete array->RemoveAt(index[i]);
3451 continue;
3452 }
e43c066c 3453 }
b9671574 3454 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3455 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
628e7bb0 3456 if (!shortbest) accepted++;
e43c066c 3457 //
3458 newarray->AddLast(array->RemoveAt(index[i]));
a5d0fae6 3459 for (Int_t j=0;j<6;j++){
3460 if (nz[j]==0){
3461 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3462 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3463 ny[j] = track->GetNy(j);
3464 nz[j] = track->GetNz(j);
e43c066c 3465 }
3466 }
3467 }
3468 else{
3469 delete array->RemoveAt(index[i]);
3470 }
3471 }
3472 array->Delete();
3473 delete fTrackHypothesys.RemoveAt(esdindex);
3474 fTrackHypothesys.AddAt(newarray,esdindex);
3475 }
3476 else{
3477 array->Delete();
3478 delete fTrackHypothesys.RemoveAt(esdindex);
3479 }
3480 }
3481 delete [] chi2;
3482 delete [] index;
3483}
44347160 3484//------------------------------------------------------------------------
15dd636f 3485AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
e43c066c 3486{
3487 //-------------------------------------------------------------
3488 // try to find best hypothesy
3489 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3490 //-------------------------------------------------------------
3491 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3492 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3493 if (!array) return 0;
3494 Int_t entries = array->GetEntriesFast();
3495 if (!entries) return 0;
3496 Float_t minchi2 = 100000;
15dd636f 3497 AliITStrackMI * besttrack=0;
e43c066c 3498 //
15dd636f 3499 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3500 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
afd25725 3501 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3502 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
e43c066c 3503 //
3504 for (Int_t i=0;i<entries;i++){
15dd636f 3505 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
e43c066c 3506 if (!track) continue;
628e7bb0 3507 Float_t sigmarfi,sigmaz;
3508 GetDCASigma(track,sigmarfi,sigmaz);
b9671574 3509 track->SetDnorm(0,sigmarfi);
3510 track->SetDnorm(1,sigmaz);
628e7bb0 3511 //
b9671574 3512 track->SetChi2MIP(1,1000000);
3513 track->SetChi2MIP(2,1000000);
3514 track->SetChi2MIP(3,1000000);
e43c066c 3515 //
3516 // backtrack
15dd636f 3517 backtrack = new(backtrack) AliITStrackMI(*track);
afd25725 3518 if (track->GetConstrain()) {
e50912db 3519 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
333d86cb 3520 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3521 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3522 }
6c94f330 3523 backtrack->ResetCovariance(10.);
628e7bb0 3524 }else{
6c94f330 3525 backtrack->ResetCovariance(10.);
628e7bb0 3526 }
e43c066c 3527 backtrack->ResetClusters();
628e7bb0 3528
e43c066c 3529 Double_t x = original->GetX();
3530 if (!RefitAt(x,backtrack,track)) continue;
628e7bb0 3531 //
b9671574 3532 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
e43c066c 3533 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
44347160 3534 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
b9671574 3535 track->SetChi22(GetMatchingChi2(backtrack,original));
e43c066c 3536
b9671574 3537 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3538 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3539 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
e43c066c 3540
3541
44347160 3542 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
e43c066c 3543 //
3544 //forward track - without constraint
15dd636f 3545 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
e43c066c 3546 forwardtrack->ResetClusters();
3547 x = track->GetX();
628e7bb0 3548 RefitAt(x,forwardtrack,track);
b9671574 3549 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
44347160 3550 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3551 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
e43c066c 3552
791f9a2a 3553 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3554 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
b9671574 3555 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3556 forwardtrack->SetD(0,track->GetD(0));
3557 forwardtrack->SetD(1,track->GetD(1));
e43c066c 3558 {
3559 Int_t list[6];
00a7cc50 3560 AliITSRecPoint* clist[6];
b9671574 3561 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3562 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
e43c066c 3563 }
3564
b9671574 3565 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
44347160 3566 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3567 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
b9671574 3568 track->SetChi2MIP(3,1000);
e43c066c 3569 continue;
3570 }
b9671574 3571 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
e43c066c 3572 //
3573 for (Int_t ichi=0;ichi<5;ichi++){
b9671574 3574 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
e43c066c 3575 }
3576 if (chi2 < minchi2){
15dd636f 3577 //besttrack = new AliITStrackMI(*forwardtrack);
99c2ee26 3578 besttrack = track;
e43c066c 3579 besttrack->SetLabel(track->GetLabel());
b9671574 3580 besttrack->SetFakeRatio(track->GetFakeRatio());
e43c066c 3581 minchi2 = chi2;
791f9a2a 3582 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3583 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
b9671574 3584 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
e43c066c 3585 }
3586 }
3587 delete backtrack;
3588 delete forwardtrack;
9cd19c63 3589
3590 if (!besttrack) return 0;
3591
e43c066c 3592 Int_t accepted=0;
3593 for (Int_t i=0;i<entries;i++){
fddf8459 3594 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3595
e43c066c 3596 if (!track) continue;
628e7bb0 3597
44347160 3598 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
e43c066c 3599 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
b9671574 3600 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3601 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
628e7bb0 3602 delete array->RemoveAt(i);
3603 continue;
3604 }
e43c066c 3605 }
3606 else{
3607 accepted++;
3608 }
3609 }
3610 //
3611 array->Compress();
3612 SortTrackHypothesys(esdindex,checkmax,1);
fddf8459 3613
e43c066c 3614 array = (TObjArray*) fTrackHypothesys.At(esdindex);
afdf9b9c 3615 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
15dd636f 3616 besttrack = (AliITStrackMI*)array->At(0);
e43c066c 3617 if (!besttrack) return 0;
b9671574 3618 besttrack->SetChi2MIP(8,0);
e43c066c 3619 fBestTrackIndex[esdindex]=0;
3620 entries = array->GetEntriesFast();
15dd636f 3621 AliITStrackMI *longtrack =0;
e43c066c 3622 minchi2 =1000;
b9671574 3623 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
6518a6c5 3624 for (Int_t itrack=entries-1;itrack>0;itrack--) {
15dd636f 3625 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
b9671574 3626 if (!track->GetConstrain()) continue;
3627 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3628 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3629 if (track->GetChi2MIP(0)>4.) continue;
3630 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
e43c066c 3631 longtrack =track;
3632 }
3633 //if (longtrack) besttrack=longtrack;
3634
3635 Int_t list[6];
00a7cc50 3636 AliITSRecPoint * clist[6];
e43c066c 3637 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
44347160 3638 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3639 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
e43c066c 3640 RegisterClusterTracks(besttrack,esdindex);
3641 }
3642 //
3643 //
3644 if (shared>0.0){
3645 Int_t nshared;
3646 Int_t overlist[6];
3647 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3648 if (sharedtrack>=0){
3649 //
3650 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3651 if (besttrack){
3652 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3653 }
3654 else return 0;
3655 }
3656 }
3657
3658 if (shared>2.5) return 0;
3659 if (shared>1.0) return besttrack;
3660 //
628e7bb0 3661 // Don't sign clusters if not gold track
3662 //
3663 if (!besttrack->IsGoldPrimary()) return besttrack;
b9671574 3664 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
628e7bb0 3665 //
e43c066c 3666 if (fConstraint[fPass]){
3667 //
3668 // sign clusters
3669 //
3670 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3671 for (Int_t i=0;i<6;i++){
b9671574 3672 Int_t index = besttrack->GetClIndex(i);
9fd412a8 3673 if (index<0) continue;
e43c066c 3674 Int_t ilayer = (index & 0xf0000000) >> 28;
b9671574 3675 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
00a7cc50 3676 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
e43c066c 3677 if (!c) continue;
3678 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3679 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3680 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
572f41f9 3681 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3682 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
e43c066c 3683 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3684
3685 Bool_t cansign = kTRUE;
3686 for (Int_t itrack=0;itrack<entries; itrack++){
15dd636f 3687 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
e43c066c 3688 if (!track) continue;
b9671574 3689 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
9fd412a8 3690 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
e43c066c 3691 cansign = kFALSE;
3692 break;
3693 }
3694 }
3695 if (cansign){
b9671574 3696 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3697 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
e43c066c 3698 if (!c->IsUsed()) c->Use();
3699 }
3700 }
3701 }
3702 return besttrack;
3703}
44347160 3704//------------------------------------------------------------------------
e43c066c 3705void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3706{
3707 //
3708 // get "best" hypothesys
628e7bb0 3709 //
e43c066c 3710
3711 Int_t nentries = itsTracks.GetEntriesFast();
3712 for (Int_t i=0;i<nentries;i++){
15dd636f 3713 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
e43c066c 3714 if (!track) continue;
3715 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3716 if (!array) continue;
3717 if (array->GetEntriesFast()<=0) continue;
3718 //
15dd636f 3719 AliITStrackMI* longtrack=0;
e43c066c 3720 Float_t minn=0;
3721 Float_t maxchi2=1000;
3722 for (Int_t j=0;j<array->GetEntriesFast();j++){
a5d0fae6 3723 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3724 if (!trackHyp) continue;
3725 if (trackHyp->GetGoldV0()) {
3726 longtrack = trackHyp; //gold V0 track taken
628e7bb0 3727 break;
3728 }
a5d0fae6 3729 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3730 Float_t chi2 = trackHyp->GetChi2MIP(0);
628e7bb0 3731 if (fAfterV0){
a5d0fae6 3732 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
628e7bb0 3733 }
a5d0fae6 3734 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
628e7bb0 3735 //
3736 if (chi2 > maxchi2) continue;
a5d0fae6 3737 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
628e7bb0 3738 maxchi2 = chi2;
a5d0fae6 3739 longtrack=trackHyp;
e43c066c 3740 }
628e7bb0 3741 //
3742 //
3743 //
15dd636f 3744 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
e43c066c 3745 if (!longtrack) {longtrack = besttrack;}
3746 else besttrack= longtrack;
628e7bb0 3747 //
afd25725 3748 if (besttrack) {
e43c066c 3749 Int_t list[6];
00a7cc50 3750 AliITSRecPoint * clist[6];
e43c066c 3751 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3752 //
b9671574 3753 track->SetNUsed(shared);
3754 track->SetNSkipped(besttrack->GetNSkipped());
3755 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
afd25725 3756 if (shared>0) {
3757 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
e43c066c 3758 Int_t nshared;
3759 Int_t overlist[6];
3760 //
3761 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3762 //if (sharedtrack==-1) sharedtrack=0;
afd25725 3763 if (sharedtrack>=0) {
e43c066c 3764 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3765 }
3766 }
afd25725 3767 if (besttrack&&fAfterV0) {
628e7bb0 3768 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3769 }
9cd19c63 3770 if (besttrack) {
3771 if (fConstraint[fPass]) UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3772 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3773 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3774 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3775 }
3776 }
e43c066c 3777 }
3778 }
3779}
44347160 3780//------------------------------------------------------------------------
15dd636f 3781void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
e43c066c 3782 //--------------------------------------------------------------------
3783 //This function "cooks" a track label. If label<0, this track is fake.
3784 //--------------------------------------------------------------------
3785 Int_t tpcLabel=-1;
3786
d3dbf323 3787 if (track->GetESDtrack()){
3788 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3789 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3790 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3791 }
b9671574 3792 track->SetChi2MIP(9,0);
e43c066c 3793 Int_t nwrong=0;
3794 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3795 Int_t cindex = track->GetClusterIndex(i);
3796 Int_t l=(cindex & 0xf0000000) >> 28;
00a7cc50 3797 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
e43c066c 3798 Int_t isWrong=1;
3799 for (Int_t ind=0;ind<3;ind++){
84a21804 3800 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
bf1cd124 3801 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
e43c066c 3802 }
b9671574 3803 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
e43c066c 3804 nwrong+=isWrong;
3805 }
7d43b4c4 3806 Int_t nclusters = track->GetNumberOfClusters();
3807 if (nclusters > 0) //PH Some tracks don't have any cluster
3808 track->SetFakeRatio(double(nwrong)/double(nclusters));
84a21804 3809 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3810 track->SetLabel(-tpcLabel);
3811 } else {
3812 track->SetLabel(tpcLabel);
e43c066c 3813 }
9be1d1c7 3814 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
e43c066c 3815
3816}
44347160 3817//------------------------------------------------------------------------
88e3c173 3818void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
3819 //
3820 // Fill the dE/dx in this track
e43c066c 3821 //
b9671574 3822 track->SetChi2MIP(9,0);
e43c066c 3823 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3824 Int_t cindex = track->GetClusterIndex(i);
3825 Int_t l=(cindex & 0xf0000000) >> 28;
00a7cc50 3826 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
b9671574 3827 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
e43c066c 3828 Int_t isWrong=1;
3829 for (Int_t ind=0;ind<3;ind++){
3830 if (cl->GetLabel(ind)==lab) isWrong=0;
3831 }
b9671574 3832 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
e43c066c 3833 }
e43c066c 3834 Double_t low=0.;
3835 Double_t up=0.51;
09cf9d66 3836 track->CookdEdx(low,up);
e43c066c 3837}
44347160 3838//------------------------------------------------------------------------
e50912db 3839void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
e43c066c 3840 //
a70ed6ad 3841 // Create some arrays
e43c066c 3842 //
afd25725 3843 if (fCoefficients) delete []fCoefficients;
3844 fCoefficients = new Float_t[ntracks*48];
3845 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
e43c066c 3846}
44347160 3847//------------------------------------------------------------------------
00a7cc50 3848Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
e43c066c 3849{
3850 //
a70ed6ad 3851 // Compute predicted chi2
e43c066c 3852 //
242bd643 3853 // Take into account the mis-alignment (bring track to cluster plane)
3854 Double_t xTrOrig=track->GetX();
3855 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
d9ead1a0 3856 Float_t erry,errz,covyz;
e43c066c 3857 Float_t theta = track->GetTgl();
3858 Float_t phi = track->GetSnp();
242bd643 3859 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
d9ead1a0 3860 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4187a601 3861 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4187a601 3862 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
d9ead1a0 3863 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
1c97ce2f 3864 // Bring the track back to detector plane in ideal geometry
3865 // [mis-alignment will be accounted for in UpdateMI()]
9be1d1c7 3866 if (!track->Propagate(xTrOrig)) return 1000.;
e43c066c 3867 Float_t ny,nz;
e50912db 3868 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
e43c066c 3869 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3870 if (delta>1){
3871 chi2+=0.5*TMath::Min(delta/2,2.);
3872 chi2+=2.*cluster->GetDeltaProbability();
3873 }
3874 //
b9671574 3875 track->SetNy(layer,ny);
3876 track->SetNz(layer,nz);
3877 track->SetSigmaY(layer,erry);
3878 track->SetSigmaZ(layer, errz);
d9ead1a0 3879 track->SetSigmaYZ(layer,covyz);
e43c066c 3880 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
60e55aee 3881 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
e43c066c 3882 return chi2;
3883
3884}
e50912db 3885//------------------------------------------------------------------------
3886Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
e43c066c 3887{
3888 //
a70ed6ad 3889 // Update ITS track
e43c066c 3890 //
3891 Int_t layer = (index & 0xf0000000) >> 28;
b9671574 3892 track->SetClIndex(layer, index);
572f41f9 3893 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3894 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3895 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3896 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3897 }
e43c066c 3898 }
a504d56f 3899
88e3c173 3900 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
a504d56f 3901
afd25725 3902
1c97ce2f 3903 // Take into account the mis-alignment (bring track to cluster plane)
3904 Double_t xTrOrig=track->GetX();
4187a601 3905 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3906 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3907 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3908 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
fddf8459 3909
9be1d1c7 3910 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
1f3e997f 3911
150f264c 3912 AliCluster c(*cl);
3913 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3914 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
d9ead1a0 3915 c.SetSigmaYZ(track->GetSigmaYZ(layer));
150f264c 3916
1c97ce2f 3917
3918 Int_t updated = track->UpdateMI(&c,chi2,index);
1c97ce2f 3919 // Bring the track back to detector plane in ideal geometry
9be1d1c7 3920 if (!track->Propagate(xTrOrig)) return 0;
d9ead1a0 3921
4187a601 3922 if(!updated) AliDebug(2,"update failed");
1c97ce2f 3923 return updated;
e43c066c 3924}
150f264c 3925
44347160 3926//------------------------------------------------------------------------
a70ed6ad 3927void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
628e7bb0 3928{
3929 //
3930 //DCA sigmas parameterization
3931 //to be paramterized using external parameters in future
3932 //
3933 //
63126a46 3934 Double_t curv=track->GetC();
3935 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
3936 sigmaz = 0.0110+4.37*TMath::Abs(curv);
628e7bb0 3937}
e50912db 3938//------------------------------------------------------------------------
a70ed6ad 3939void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
e43c066c 3940{
3941 //
a70ed6ad 3942 // Clusters from delta electrons?
e43c066c 3943 //
a70ed6ad 3944 Int_t entries = clusterArray->GetEntriesFast();
e43c066c 3945 if (entries<4) return;
a70ed6ad 3946 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
e43c066c 3947 Int_t layer = cluster->GetLayer();
3948 if (layer>1) return;
3949 Int_t index[10000];
3950 Int_t ncandidates=0;
3951 Float_t r = (layer>0)? 7:4;
3952 //
3953 for (Int_t i=0;i<entries;i++){
a70ed6ad 3954 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
e43c066c 3955 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3956 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3957 index[ncandidates] = i; //candidate to belong to delta electron track
3958 ncandidates++;
3959 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3960 cl0->SetDeltaProbability(1);
3961 }
3962 }
3963 //
3964 //
3965 //
3966 for (Int_t i=0;i<ncandidates;i++){
a70ed6ad 3967 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
e43c066c 3968 if (cl0->GetDeltaProbability()>0.8) continue;
3969 //
3970 Int_t ncl = 0;
3971 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3972 sumy=sumz=sumy2=sumyz=sumw=0.0;
3973 for (Int_t j=0;j<ncandidates;j++){
3974 if (i==j) continue;
a70ed6ad 3975 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
e43c066c 3976 //
3977 Float_t dz = cl0->GetZ()-cl1->GetZ();
3978 Float_t dy = cl0->GetY()-cl1->GetY();
3979 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3980 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3981 y[ncl] = cl1->GetY();
3982 z[ncl] = cl1->GetZ();
3983 sumy+= y[ncl]*weight;
3984 sumz+= z[ncl]*weight;
3985 sumy2+=y[ncl]*y[ncl]*weight;
3986 sumyz+=y[ncl]*z[ncl]*weight;
3987 sumw+=weight;
3988 ncl++;
3989 }
3990 }
3991 if (ncl<4) continue;
3992 Float_t det = sumw*sumy2 - sumy*sumy;
3993 Float_t delta=1000;
3994 if (TMath::Abs(det)>0.01){
3995 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3996 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3997 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3998 }
3999 else{
4000 Float_t z0 = sumyz/sumy;
4001 delta = TMath::Abs(cl0->GetZ()-z0);
4002 }
4003 if ( delta<0.05) {
4004 cl0->SetDeltaProbability(1-20.*delta);
4005 }
4006 }
4007}
44347160 4008//------------------------------------------------------------------------
15dd636f 4009void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
e43c066c 4010{
4011 //
a70ed6ad 4012 // Update ESD track
e43c066c 4013 //
4014 track->UpdateESDtrack(flags);
b9671574 4015 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
99c2ee26 4016 if (oldtrack) delete oldtrack;
b9671574 4017 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
5668b5b2 4018 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4019 // printf("Problem\n");
4020 // }
e43c066c 4021}
44347160 4022//------------------------------------------------------------------------
81e97e0d 4023Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4024 //
ae00569a 4025 // Get nearest upper layer close to the point xr.
81e97e0d 4026 // rough approximation
4027 //
bf6adc12 4028 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
81e97e0d 4029 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4030 Int_t res =6;
4031 for (Int_t i=0;i<6;i++){
bf6adc12 4032 if (radius<kRadiuses[i]){
81e97e0d 4033 res =i;
4034 break;
4035 }
4036 }
4037 return res;
4038}
44347160 4039//------------------------------------------------------------------------
e50912db 4040void AliITStrackerMI::BuildMaterialLUT(TString material) {
4041 //--------------------------------------------------------------------
4042 // Fill a look-up table with mean material
4043 //--------------------------------------------------------------------
4044
4045 Int_t n=1000;
4046 Double_t mparam[7];
4047 Double_t point1[3],point2[3];
4048 Double_t phi,cosphi,sinphi,z;
4049 // 0-5 layers, 6 pipe, 7-8 shields
1ec3c859 4050 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4051 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
e50912db 4052
4053 Int_t ifirst=0,ilast=0;
4054 if(material.Contains("Pipe")) {
4055 ifirst=6; ilast=6;
4056 } else if(material.Contains("Shields")) {
4057 ifirst=7; ilast=8;
4058 } else if(material.Contains("Layers")) {
4059 ifirst=0; ilast=5;
4060 } else {
4061 Error("BuildMaterialLUT","Wrong layer name\n");
4062 }
4063
4064 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4065 Double_t param[5]={0.,0.,0.,0.,0.};
4066 for (Int_t i=0; i<n; i++) {
4067 phi = 2.*TMath::Pi()*gRandom->Rndm();
4068 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4069 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4070 point1[0] = rmin[imat]*cosphi;
4071 point1[1] = rmin[imat]*sinphi;
4072 point1[2] = z;
4073 point2[0] = rmax[imat]*cosphi;
4074 point2[1] = rmax[imat]*sinphi;
4075 point2[2] = z;
4076 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4077 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4078 }
4079 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4080 if(imat<=5) {
4081 fxOverX0Layer[imat] = param[1];
4082 fxTimesRhoLayer[imat] = param[0]*param[4];
4083 } else if(imat==6) {
4084 fxOverX0Pipe = param[1];
4085 fxTimesRhoPipe = param[0]*param[4];
4086 } else if(imat==7) {
4087 fxOverX0Shield[0] = param[1];
4088 fxTimesRhoShield[0] = param[0]*param[4];
4089 } else if(imat==8) {
4090 fxOverX0Shield[1] = param[1];
4091 fxTimesRhoShield[1] = param[0]*param[4];
4092 }
4093 }
4094 /*
4095 printf("%s\n",material.Data());
4096 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4097 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4098 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4099 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4100 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4101 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4102 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4103 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4104 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4105 */
4106 return;
4107}
4108//------------------------------------------------------------------------
4109Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4110 TString direction) {
4111 //-------------------------------------------------------------------
4112 // Propagate beyond beam pipe and correct for material
4113 // (material budget in different ways according to fUseTGeo value)
7131593b 4114 // Add time if going outward (PropagateTo or PropagateToTGeo)
e50912db 4115 //-------------------------------------------------------------------
4116
4117 // Define budget mode:
4118 // 0: material from AliITSRecoParam (hard coded)
9be1d1c7 4119 // 1: material from TGeo in one step (on the fly)
e50912db 4120 // 2: material from lut
9be1d1c7 4121 // 3: material from TGeo in one step (same for all hypotheses)
e50912db 4122 Int_t mode;
4123 switch(fUseTGeo) {
4124 case 0:
4125 mode=0;
4126 break;
4127 case 1:
4128 mode=1;
4129 break;
4130 case 2:
4131 mode=2;
4132 break;
4133 case 3:
4134 if(fTrackingPhase.Contains("Clusters2Tracks"))
4135 { mode=3; } else { mode=1; }
4136 break;
4137 case 4:
4138 if(fTrackingPhase.Contains("Clusters2Tracks"))
4139 { mode=3; } else { mode=2; }
4140 break;
4141 default:
4142 mode=0;
4143 break;
4144 }
4145 if(fTrackingPhase.Contains("Default")) mode=0;
4146
4147 Int_t index=fCurrentEsdTrack;
4148
4149 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4150 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
70e83b26 4151 Double_t xToGo;
4152 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
e50912db 4153
4154 Double_t xOverX0,x0,lengthTimesMeanDensity;
e50912db 4155
4156 switch(mode) {
4157 case 0:
4158 xOverX0 = AliITSRecoParam::GetdPipe();
4159 x0 = AliITSRecoParam::GetX0Be();
4160 lengthTimesMeanDensity = xOverX0*x0;
9be1d1c7 4161 lengthTimesMeanDensity *= dir;
7131593b 4162 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4163 break;
4164 case 1:
4165 if (!t->PropagateToTGeo(xToGo,1)) return 0;
e50912db 4166 break;
4167 case 2:
4168 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4169 xOverX0 = fxOverX0Pipe;
4170 lengthTimesMeanDensity = fxTimesRhoPipe;
9be1d1c7 4171 lengthTimesMeanDensity *= dir;
7131593b 4172 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4173 break;
4174 case 3:
4175 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4176 if(fxOverX0PipeTrks[index]<0) {
4177 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4178 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
60e55aee 4179 ((1.-t->GetSnp())*(1.+t->GetSnp())));
e50912db 4180 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4181 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4182 return 1;
4183 }
4184 xOverX0 = fxOverX0PipeTrks[index];
4185 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
9be1d1c7 4186 lengthTimesMeanDensity *= dir;
7131593b 4187 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4188 break;
4189 }
4190
e50912db 4191 return 1;
4192}
4193//------------------------------------------------------------------------
4194Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4195 TString shield,
4196 TString direction) {
4197 //-------------------------------------------------------------------
4198 // Propagate beyond SPD or SDD shield and correct for material
4199 // (material budget in different ways according to fUseTGeo value)
7131593b 4200 // Add time if going outward (PropagateTo or PropagateToTGeo)
e50912db 4201 //-------------------------------------------------------------------
4202
4203 // Define budget mode:
4204 // 0: material from AliITSRecoParam (hard coded)
9be1d1c7 4205 // 1: material from TGeo in steps of X cm (on the fly)
4206 // X = AliITSRecoParam::GetStepSizeTGeo()
e50912db 4207 // 2: material from lut
9be1d1c7 4208 // 3: material from TGeo in one step (same for all hypotheses)
e50912db 4209 Int_t mode;
4210 switch(fUseTGeo) {
4211 case 0:
4212 mode=0;
4213 break;
4214 case 1:
4215 mode=1;
4216 break;
4217 case 2:
4218 mode=2;
4219 break;
4220 case 3:
4221 if(fTrackingPhase.Contains("Clusters2Tracks"))
4222 { mode=3; } else { mode=1; }
4223 break;
4224 case 4:
4225 if(fTrackingPhase.Contains("Clusters2Tracks"))
4226 { mode=3; } else { mode=2; }
4227 break;
4228 default:
4229 mode=0;
4230 break;
4231 }
4232 if(fTrackingPhase.Contains("Default")) mode=0;
4233
4234 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4235 Double_t rToGo;
4236 Int_t shieldindex=0;
4237 if (shield.Contains("SDD")) { // SDDouter
4238 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4239 shieldindex=1;
4240 } else if (shield.Contains("SPD")) { // SPDouter
4241 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4242 shieldindex=0;
4243 } else {
4244 Error("CorrectForShieldMaterial"," Wrong shield name\n");
e50912db 4245 return 0;
4246 }
67178f34 4247
4248 // do nothing if we are already beyond the shield
4249 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4250 if(dir<0 && rTrack > rToGo) return 1; // going outward
4251 if(dir>0 && rTrack < rToGo) return 1; // going inward
4252
4253
70e83b26 4254 Double_t xToGo;
4255 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
e50912db 4256
4257 Int_t index=2*fCurrentEsdTrack+shieldindex;
4258
4259 Double_t xOverX0,x0,lengthTimesMeanDensity;
9be1d1c7 4260 Int_t nsteps=1;
e50912db 4261
4262 switch(mode) {
4263 case 0:
4264 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4265 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4266 lengthTimesMeanDensity = xOverX0*x0;
9be1d1c7 4267 lengthTimesMeanDensity *= dir;
7131593b 4268 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4269 break;
4270 case 1:
9be1d1c7 4271 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4272 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
e50912db 4273 break;
4274 case 2:
4275 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4276 xOverX0 = fxOverX0Shield[shieldindex];
4277 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
9be1d1c7 4278 lengthTimesMeanDensity *= dir;
7131593b 4279 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4280 break;
4281 case 3:
4282 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4283 if(fxOverX0ShieldTrks[index]<0) {
4284 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4285 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
60e55aee 4286 ((1.-t->GetSnp())*(1.+t->GetSnp())));
e50912db 4287 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4288 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4289 return 1;
4290 }
4291 xOverX0 = fxOverX0ShieldTrks[index];
4292 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
9be1d1c7 4293 lengthTimesMeanDensity *= dir;
7131593b 4294 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4295 break;
4296 }
4297
e50912db 4298 return 1;
4299}
4300//------------------------------------------------------------------------
4301Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4302 Int_t layerindex,
4303 Double_t oldGlobXYZ[3],
4304 TString direction) {
4305 //-------------------------------------------------------------------
4306 // Propagate beyond layer and correct for material
4307 // (material budget in different ways according to fUseTGeo value)
7131593b 4308 // Add time if going outward (PropagateTo or PropagateToTGeo)
e50912db 4309 //-------------------------------------------------------------------
4310
4311 // Define budget mode:
4312 // 0: material from AliITSRecoParam (hard coded)
9be1d1c7 4313 // 1: material from TGeo in stepsof X cm (on the fly)
4314 // X = AliITSRecoParam::GetStepSizeTGeo()
e50912db 4315 // 2: material from lut
9be1d1c7 4316 // 3: material from TGeo in one step (same for all hypotheses)
e50912db 4317 Int_t mode;
4318 switch(fUseTGeo) {
4319 case 0:
4320 mode=0;
4321 break;
4322 case 1:
4323 mode=1;
4324 break;
4325 case 2:
4326 mode=2;
4327 break;
4328 case 3:
4329 if(fTrackingPhase.Contains("Clusters2Tracks"))
4330 { mode=3; } else { mode=1; }
4331 break;
4332 case 4:
4333 if(fTrackingPhase.Contains("Clusters2Tracks"))
4334 { mode=3; } else { mode=2; }
4335 break;
4336 default:
4337 mode=0;
4338 break;
4339 }
4340 if(fTrackingPhase.Contains("Default")) mode=0;
4341
4342 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4343
4344 Double_t r=fgLayers[layerindex].GetR();
4345 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4346
4347 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
70e83b26 4348 Double_t xToGo;
4349 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
e50912db 4350
4351 Int_t index=6*fCurrentEsdTrack+layerindex;
4352
e50912db 4353
4354 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
9be1d1c7 4355 Int_t nsteps=1;
7131593b 4356
4357 // back before material (no correction)
d634c32b 4358 Double_t rOld,xOld;
7131593b 4359 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4360 if (!t->GetLocalXat(rOld,xOld)) return 0;
4361 if (!t->Propagate(xOld)) return 0;
e50912db 4362
4363 switch(mode) {
4364 case 0:
4365 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4366 lengthTimesMeanDensity = xOverX0*x0;
9be1d1c7 4367 lengthTimesMeanDensity *= dir;
7131593b 4368 // Bring the track beyond the material
4369 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4370 break;
4371 case 1:
9be1d1c7 4372 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
7131593b 4373 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
e50912db 4374 break;
4375 case 2:
4376 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4377 xOverX0 = fxOverX0Layer[layerindex];
4378 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
9be1d1c7 4379 lengthTimesMeanDensity *= dir;
7131593b 4380 // Bring the track beyond the material
4381 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4382 break;
4383 case 3:
4384 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
7131593b 4385 if(fxOverX0LayerTrks[index]<0) {
4386 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4387 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
e50912db 4388 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
60e55aee 4389 ((1.-t->GetSnp())*(1.+t->GetSnp())));
7131593b 4390 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4391 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4392 return 1;
e50912db 4393 }
4394 xOverX0 = fxOverX0LayerTrks[index];
4395 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
9be1d1c7 4396 lengthTimesMeanDensity *= dir;
7131593b 4397 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
e50912db 4398 break;
4399 }
4400
7131593b 4401
e50912db 4402 return 1;
4403}
4404//------------------------------------------------------------------------
4405void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4406 //-----------------------------------------------------------------
4407 // Initialize LUT for storing material for each prolonged track
4408 //-----------------------------------------------------------------
4409 fxOverX0PipeTrks = new Float_t[ntracks];
4410 fxTimesRhoPipeTrks = new Float_t[ntracks];
4411 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4412 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4413 fxOverX0LayerTrks = new Float_t[ntracks*6];
4414 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4415
4416 for(Int_t i=0; i<ntracks; i++) {
4417 fxOverX0PipeTrks[i] = -1.;
4418 fxTimesRhoPipeTrks[i] = -1.;
4419 }
4420 for(Int_t j=0; j<ntracks*2; j++) {
4421 fxOverX0ShieldTrks[j] = -1.;
4422 fxTimesRhoShieldTrks[j] = -1.;
4423 }
4424 for(Int_t k=0; k<ntracks*6; k++) {
4425 fxOverX0LayerTrks[k] = -1.;
4426 fxTimesRhoLayerTrks[k] = -1.;
4427 }
4428
4429 fNtracks = ntracks;
4430
4431 return;
4432}
4433//------------------------------------------------------------------------
4434void AliITStrackerMI::DeleteTrksMaterialLUT() {
4435 //-----------------------------------------------------------------
4436 // Delete LUT for storing material for each prolonged track
4437 //-----------------------------------------------------------------
4438 if(fxOverX0PipeTrks) {
4439 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4440 }
4441 if(fxOverX0ShieldTrks) {
4442 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4443 }
4444
4445 if(fxOverX0LayerTrks) {
4446 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4447 }
4448 if(fxTimesRhoPipeTrks) {
4449 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4450 }
4451 if(fxTimesRhoShieldTrks) {
4452 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4453 }
4454 if(fxTimesRhoLayerTrks) {
4455 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4456 }
4457 return;
4458}
4459//------------------------------------------------------------------------
4fa7f7d1 4460void AliITStrackerMI::SetForceSkippingOfLayer() {
4461 //-----------------------------------------------------------------
4462 // Check if we are forced to skip layers
4463 // either we set to skip them in RecoParam
4464 // or they were off during data-taking
4465 //-----------------------------------------------------------------
4466
4467 const AliEventInfo *eventInfo = GetEventInfo();
6262dd3d 4468
4fa7f7d1 4469 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4470 fForceSkippingOfLayer[l] = 0;
4471 // check reco param
4472 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4473 // check run info
4474
87b4605f 4475 if(eventInfo &&
4476 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4fa7f7d1 4477 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4478 if(l==0 || l==1) {
4479 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4480 } else if(l==2 || l==3) {
4481 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4482 } else {
4483 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4484 }
4485 }
4486 }
4487 return;
4488}
4489//------------------------------------------------------------------------
a70ed6ad 4490Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
6518a6c5 4491 Int_t ilayer,Int_t idet) const {
2755f080 4492 //-----------------------------------------------------------------
4493 // This method is used to decide whether to allow a prolongation
4494 // without clusters, because we want to skip the layer.
4495 // In this case the return value is > 0:
4496 // return 1: the user requested to skip a layer
371164ec 4497 // return 2: track outside z acceptance
2755f080 4498 //-----------------------------------------------------------------
e43c066c 4499
4fa7f7d1 4500 if (ForceSkippingOfLayer(ilayer)) return 1;
51ad6848 4501
371164ec 4502 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4503
4504 if (idet<0 && // out in z
4505 ilayer>innerLayCanSkip &&
4506 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
2755f080 4507 // check if track will cross SPD outer layer
4508 Double_t phiAtSPD2,zAtSPD2;
6518a6c5 4509 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
2755f080 4510 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4511 }
371164ec 4512 return 2; // always allow skipping, changed on 05.11.2009
2755f080 4513 }
51ad6848 4514
2755f080 4515 return 0;
4516}
4517//------------------------------------------------------------------------
23197852 4518Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
ae00569a 4519 Int_t ilayer,Int_t idet,
1c97ce2f 4520 Double_t dz,Double_t dy,
23197852 4521 Bool_t noClusters) const {
6518a6c5 4522 //-----------------------------------------------------------------
4523 // This method is used to decide whether to allow a prolongation
ae00569a 4524 // without clusters, because there is a dead zone in the road.
6518a6c5 4525 // In this case the return value is > 0:
4526 // return 1: dead zone at z=0,+-7cm in SPD
bf1cd124 4527 // This method assumes that fSPDdetzcentre is ordered from -z to +z
23197852 4528 // return 2: all road is "bad" (dead or noisy) from the OCDB
f8720bda 4529 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4530 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
6518a6c5 4531 //-----------------------------------------------------------------
4532
4533 // check dead zones at z=0,+-7cm in the SPD
4534 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4535 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4536 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4537 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4538 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4539 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4540 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4541 for (Int_t i=0; i<3; i++)
4187a601 4542 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
bf1cd124 4543 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4544 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4187a601 4545 }
6518a6c5 4546 }
4547
23197852 4548 // check bad zones from OCDB
4549 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
ae00569a 4550
23197852 4551 if (idet<0) return 0;
ae00569a 4552
23197852 4553 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4554
23197852 4555 Int_t detType=-1;
4556 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4557 if (ilayer==0 || ilayer==1) { // ---------- SPD
4558 detType = 0;
4559 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4560 detType = 1;
4561 detSizeFactorX *= 2.;
4562 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4563 detType = 2;
4564 }
a70ed6ad 4565 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
23197852 4566 if (detType==2) segm->SetLayer(ilayer+1);
4567 Float_t detSizeX = detSizeFactorX*segm->Dx();
4568 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4569
4570 // check if the road overlaps with bad chips
4571 Float_t xloc,zloc;
e9b15b0c 4572 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
1c97ce2f 4573 Float_t zlocmin = zloc-dz;
4574 Float_t zlocmax = zloc+dz;
4575 Float_t xlocmin = xloc-dy;
4576 Float_t xlocmax = xloc+dy;
23197852 4577 Int_t chipsInRoad[100];
4578
4187a601 4579 // check if road goes out of detector
4580 Bool_t touchNeighbourDet=kFALSE;
6262dd3d 4581 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4582 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4583 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4584 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4187a601 4585 AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
4586
4587 // check if this detector is bad
4588 if (det.IsBad()) {
4589 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4590 if(!touchNeighbourDet) {
4591 return 2; // all detectors in road are bad
4592 } else {
4593 return 3; // at least one is bad
4594 }
4595 }
4596
23197852 4597 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4187a601 4598 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
23197852 4599 if (!nChipsInRoad) return 0;
4600
4601 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4602 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4603 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4187a601 4604 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
23197852 4605 if (det.IsChipBad(chipsInRoad[iCh])) {
4606 anyBad=kTRUE;
4607 } else {
4608 anyGood=kTRUE;
4609 }
6518a6c5 4610 }
4611
4187a601 4612 if (!anyGood) {
4613 if(!touchNeighbourDet) {
4614 AliDebug(2,"all bad in road");
4615 return 2; // all chips in road are bad
4616 } else {
4617 return 3; // at least a bad chip in road
4618 }
4619 }
23197852 4620
4187a601 4621 if (anyBad) {
4622 AliDebug(2,"at least a bad in road");
4623 return 3; // at least a bad chip in road
4624 }
23197852 4625
4626
4627 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
23197852 4628 || !noClusters) return 0;
4629
4630 // There are no clusters in road: check if there is at least
6a0fdcbd 4631 // a bad SPD pixel or SDD anode or SSD strips on both sides
23197852 4632
4187a601 4633 Int_t idetInITS=idet;
4634 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
23197852 4635
4187a601 4636 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4637 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
f8720bda 4638 return 4;
4187a601 4639 }
4640 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
23197852 4641
6518a6c5 4642 return 0;
4643}
4644//------------------------------------------------------------------------
ae00569a 4645Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
a70ed6ad 4646 const AliITStrackMI *track,
ae00569a 4647 Float_t &xloc,Float_t &zloc) const {
4648 //-----------------------------------------------------------------
4649 // Gives position of track in local module ref. frame
4650 //-----------------------------------------------------------------
4651
4652 xloc=0.;
4653 zloc=0.;
4654
371164ec 4655 if(idet<0) return kTRUE; // track out of z acceptance of layer
ae00569a 4656
4657 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4658
4659 Int_t lad = Int_t(idet/ndet) + 1;
4a66240a 4660
ae00569a 4661 Int_t det = idet - (lad-1)*ndet + 1;
4662
4663 Double_t xyzGlob[3],xyzLoc[3];
4664
1c97ce2f 4665 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4666 // take into account the misalignment: xyz at real detector plane
d1181f10 4667 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
ae00569a 4668
d1181f10 4669 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
ae00569a 4670
4671 xloc = (Float_t)xyzLoc[0];
4672 zloc = (Float_t)xyzLoc[2];
4673
4674 return kTRUE;
4675}
4676//------------------------------------------------------------------------
879cdb02 4677//------------------------------------------------------------------------
a70ed6ad 4678Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
0e91092f 4679//
4680// Method to be optimized further:
4681// Aim: decide whether a track can be used for PlaneEff evaluation
4682// the decision is taken based on the track quality at the layer under study
4683// no information on the clusters on this layer has to be used
4684// The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4685// the cut is done on number of sigmas from the boundaries
4686//
4687// Input: Actual track, layer [0,5] under study
4688// Output: none
4689// Return: kTRUE if this is a good track
7167ae53 4690//
4691// it will apply a pre-selection to obtain good quality tracks.
4692// Here also you will have the possibility to put a control on the
4693// impact point of the track on the basic block, in order to exclude border regions
4694// this will be done by calling a proper method of the AliITSPlaneEff class.
4695//
4696// input: AliITStrackMI* track, ilayer= layer number [0,5]
0e91092f 4697// return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
0ed58a47 4698//
4699 Int_t index[AliITSgeomTGeo::kNLayers];
4700 Int_t k;
4701 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4702 //
4703 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4704 index[k]=clusters[k];
4705 }
4706
0e91092f 4707 if(!fPlaneEff)
a31a139f 4708 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
8a32c4b5 4709 AliITSlayer &layer=fgLayers[ilayer];
4710 Double_t r=layer.GetR();
8a32c4b5 4711 AliITStrackMI tmp(*track);
4712
0ed58a47 4713// require a minimal number of cluster in other layers and eventually clusters in closest layers
061c42a0 4714 Int_t ncl_out=0; Int_t ncl_in=0;
4715 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
879cdb02 4716 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
0ed58a47 4717 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
879cdb02 4718 // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4719if(index[lay]>=0)ncl_out++;
0ed58a47 4720 }
061c42a0 4721 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4722 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4723 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
879cdb02 4724 if (index[lay]>=0) ncl_in++;
061c42a0 4725 }
63cce7b6 4726 Int_t ncl=ncl_out+ncl_in;
0ed58a47 4727 Bool_t nextout = kFALSE;
4728 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
9fd412a8 4729 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
0ed58a47 4730 Bool_t nextin = kFALSE;
4731 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4732 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
061c42a0 4733 // maximum number of missing clusters allowed in outermost layers
4734 if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4735 return kFALSE;
4736 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4737 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
0ed58a47 4738 return kFALSE;
4739 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4740 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4741 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4742 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4743
8a32c4b5 4744// detector number
4745 Double_t phi,z;
4746 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4747 Int_t idet=layer.FindDetectorIndex(phi,z);
4748 if(idet<0) { AliInfo(Form("cannot find detector"));
4749 return kFALSE;}
4750
4751 // here check if it has good Chi Square.
4752
4753 //propagate to the intersection with the detector plane
4754 const AliITSdetector &det=layer.GetDetector(idet);
4755 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4756
4757 Float_t locx; //
4758 Float_t locz; //
0e91092f 4759 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
8a32c4b5 4760 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4761 if(key>fPlaneEff->Nblock()) return kFALSE;
4762 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4763 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
8a32c4b5 4764 //***************
0e91092f 4765 // DEFINITION OF SEARCH ROAD FOR accepting a track
8a32c4b5 4766 //
061c42a0 4767 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4768 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
0e91092f 4769 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4770 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
879cdb02 4771 // done for RecPoints
8a32c4b5 4772
1c97ce2f 4773 // exclude tracks at boundary between detectors
0e91092f 4774 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
8a32c4b5 4775 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
a31a139f 4776 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
0e91092f 4777 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4778 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
0e91092f 4779 if ( (locx-dx < blockXmn+boundaryWidth) ||
4780 (locx+dx > blockXmx-boundaryWidth) ||
4781 (locz-dz < blockZmn+boundaryWidth) ||
4782 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
8a32c4b5 4783 return kTRUE;
7167ae53 4784}
4785//------------------------------------------------------------------------
a70ed6ad 4786void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
7167ae53 4787//
4788// This Method has to be optimized! For the time-being it uses the same criteria
4789// as those used in the search of extra clusters for overlapping modules.
4790//
4791// Method Purpose: estabilish whether a track has produced a recpoint or not
4792// in the layer under study (For Plane efficiency)
4793//
4794// inputs: AliITStrackMI* track (pointer to a usable track)
4795// outputs: none
4796// side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4797// traversed by this very track. In details:
4798// - if a cluster can be associated to the track then call
4799// AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4800// - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4801//
a31a139f 4802 if(!fPlaneEff)
4803 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
7167ae53 4804 AliITSlayer &layer=fgLayers[ilayer];
4805 Double_t r=layer.GetR();
7167ae53 4806 AliITStrackMI tmp(*track);
4807
4808// detector number
4809 Double_t phi,z;
4810 if (!tmp.GetPhiZat(r,phi,z)) return;
4811 Int_t idet=layer.FindDetectorIndex(phi,z);
4812
8a32c4b5 4813 if(idet<0) { AliInfo(Form("cannot find detector"));
4814 return;}
7167ae53 4815
7167ae53 4816
0e91092f 4817//propagate to the intersection with the detector plane
7167ae53 4818 const AliITSdetector &det=layer.GetDetector(idet);
4819 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
8a32c4b5 4820
1c97ce2f 4821
0e91092f 4822//***************
4823// DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4824//
4825 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4826 TMath::Sqrt(tmp.GetSigmaZ2() +
4827 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4828 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4829 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4830 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4831 TMath::Sqrt(tmp.GetSigmaY2() +
4832 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4833 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4834 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4835
4836// road in global (rphi,z) [i.e. in tracking ref. system]
4837 Double_t zmin = tmp.GetZ() - dz;
4838 Double_t zmax = tmp.GetZ() + dz;
4839 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4840 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4841
4842// select clusters in road
7167ae53 4843 layer.SelectClusters(zmin,zmax,ymin,ymax);
4844
0e91092f 4845// Define criteria for track-cluster association
7167ae53 4846 Double_t msz = tmp.GetSigmaZ2() +
4847 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4848 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4849 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4850 Double_t msy = tmp.GetSigmaY2() +
4851 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4852 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4853 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
8a32c4b5 4854 if (tmp.GetConstrain()) {
7167ae53 4855 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4856 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
8a32c4b5 4857 } else {
4858 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4859 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4860 }
7167ae53 4861 msz = 1./msz; // 1/RoadZ^2
4862 msy = 1./msy; // 1/RoadY^2
0e91092f 4863//
7167ae53 4864
4865 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4866 Int_t idetc=-1;
8a32c4b5 4867 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4868 //Double_t tolerance=0.2;
4869 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
7167ae53 4870 idetc = cl->GetDetectorIndex();
4871 if(idet!=idetc) continue;
4872 //Int_t ilay = cl->GetLayer();
4873
4874 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4875 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4876
4877 Double_t chi2=tmp.GetPredictedChi2(cl);
8a32c4b5 4878 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4879 }*/
4880 Float_t locx; //
4881 Float_t locz; //
0e91092f 4882 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
7167ae53 4883//
4884 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4885 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4886 if(key>fPlaneEff->Nblock()) return;
4887 Bool_t found=kFALSE;
8a32c4b5 4888 //if (ci>=0) {
1cc5cedc 4889 Double_t chi2;
8a32c4b5 4890 while ((cl=layer.GetNextCluster(clidx))!=0) {
4891 idetc = cl->GetDetectorIndex();
4892 if(idet!=idetc) continue;
4893 // here real control to see whether the cluster can be associated to the track.
4894 // cluster not associated to track
4895 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4896 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
4897 // calculate track-clusters chi2
1cc5cedc 4898 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4899 // in particular, the error associated to the cluster
8a32c4b5 4900 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4901 // chi2 cut
4902 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
7167ae53 4903 found=kTRUE;
8a32c4b5 4904 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
7167ae53 4905 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4906 // track->SetExtraModule(ilayer,idetExtra);
4907 }
4908 if(!fPlaneEff->UpDatePlaneEff(found,key))
8a32c4b5 4909 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5fbd4fd6 4910 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
1cc5cedc 4911 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
4912 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4913 Int_t cltype[2]={-999,-999};
879cdb02 4914 // and the module
4915
4916Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
0e91092f 4917
4918 tr[0]=locx;
4919 tr[1]=locz;
4920 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4921 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4922
3ebe30ad 4923 if (found){
4924 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4925 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4926 cltype[0]=layer.GetCluster(ci)->GetNy();
4927 cltype[1]=layer.GetCluster(ci)->GetNz();
1cc5cedc 4928
4929 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4930 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
4931 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4932 // It is computed properly by calling the method
4933 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4934 // T
4935 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4936 //if (tmp.PropagateTo(x,0.,0.)) {
4937 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4938 AliCluster c(*layer.GetCluster(ci));
4939 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4940 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
1cc5cedc 4941 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
0e91092f 4942 clu[2]=TMath::Sqrt(c.GetSigmaY2());
4943 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
1cc5cedc 4944 //}
3ebe30ad 4945 }
879cdb02 4946 // Compute the angles between the track and the module
4947 // compute the angle "in phi direction", i.e. the angle in the transverse plane
4948 // between the normal to the module and the projection (in the transverse plane) of the
4949 // track trajectory
4950 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
4951 Float_t tgl = tmp.GetTgl();
4952 Float_t phitr = tmp.GetSnp();
4953 phitr = TMath::ASin(phitr);
4954 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
4955
4956 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
4957 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
4958 Double_t alpha =0.;
4959 alpha = tmp.GetAlpha();
4960 Double_t phiglob = alpha+phitr;
4961 Double_t p[3];
4962 p[0] = TMath::Cos(phiglob);
4963 p[1] = TMath::Sin(phiglob);
4964 p[2] = tgl;
4965 TVector3 pvec(p[0],p[1],p[2]);
4966 TVector3 normvec(rot[1],rot[4],rot[7]);
4967 Double_t angle = pvec.Angle(normvec);
4968
4969 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
4970 angle *= 180./TMath::Pi();
4971
4972 //Trasverse Plane
4973 TVector3 pt(p[0],p[1],0);
4974 TVector3 normt(rot[1],rot[4],0);
4975 Double_t anglet = pt.Angle(normt);
4976
4977 Double_t phiPt = TMath::ATan2(p[1],p[0]);
4978 if(phiPt<0)phiPt+=2.*TMath::Pi();
4979 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
4980 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
4981 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
4982 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
4983 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
4984 anglet *= 180./TMath::Pi();
4985
4986 AngleModTrack[2]=(Float_t) angle;
4987 AngleModTrack[0]=(Float_t) anglet;
4988 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
4989 AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
4990 AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
4991 AngleModTrack[1]*=180./TMath::Pi(); // in degree
4992
4993 fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
5fbd4fd6 4994 }
7167ae53 4995return;
4996}