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