1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------------------------
18 // Implementation of the ITS tracker class
19 // It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
20 // and fills with them the ESD
21 // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
22 // Current support and development:
23 // Andrea Dainese, andrea.dainese@lnl.infn.it
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
27 //-------------------------------------------------------------------------
31 #include <TDatabasePDG.h>
34 #include <TTreeStream.h>
39 #include "AliGeomManager.h"
40 #include "AliITSPlaneEff.h"
41 #include "AliTrackPointArray.h"
42 #include "AliESDEvent.h"
43 #include "AliESDV0Params.h"
44 #include "AliESDtrack.h"
46 #include "AliITSChannelStatus.h"
47 #include "AliITSDetTypeRec.h"
48 #include "AliITSRecPoint.h"
49 #include "AliITSRecPointContainer.h"
50 #include "AliITSgeomTGeo.h"
51 #include "AliITSReconstructor.h"
52 #include "AliITSClusterParam.h"
53 #include "AliITSsegmentation.h"
54 #include "AliITSCalibration.h"
55 #include "AliITSPlaneEffSPD.h"
56 #include "AliITSPlaneEffSDD.h"
57 #include "AliITSPlaneEffSSD.h"
58 #include "AliITSV0Finder.h"
59 #include "AliITStrackerMI.h"
60 #include "AliMathBase.h"
64 ClassImp(AliITStrackerMI)
66 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
68 AliITStrackerMI::AliITStrackerMI():AliTracker(),
78 fLastLayerToTrackTo(0),
81 fTrackingPhase("Default"),
85 fSelectBestMIP03(kFALSE),
86 fUseImproveKalman(kFALSE),
90 fxTimesRhoPipeTrks(0),
91 fxOverX0ShieldTrks(0),
92 fxTimesRhoShieldTrks(0),
94 fxTimesRhoLayerTrks(0),
99 fSPDChipIntPlaneEff(0),
103 //Default constructor
105 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
107 fxOverX0Shield[i]=-1.;
108 fxTimesRhoShield[i]=-1.;
111 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
112 fOriginal.SetOwner();
113 for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
114 for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
115 fITSPid=new AliITSPIDResponse();
118 //------------------------------------------------------------------------
119 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
120 fI(AliITSgeomTGeo::GetNLayers()),
129 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
132 fTrackingPhase("Default"),
136 fSelectBestMIP03(kFALSE),
137 fUseImproveKalman(kFALSE),
141 fxTimesRhoPipeTrks(0),
142 fxOverX0ShieldTrks(0),
143 fxTimesRhoShieldTrks(0),
144 fxOverX0LayerTrks(0),
145 fxTimesRhoLayerTrks(0),
147 fITSChannelStatus(0),
150 fSPDChipIntPlaneEff(0),
152 //--------------------------------------------------------------------
153 //This is the AliITStrackerMI constructor
154 //--------------------------------------------------------------------
156 AliWarning("\"geom\" is actually a dummy argument !");
159 fOriginal.SetOwner();
163 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
164 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
165 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
167 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
168 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
169 Double_t poff=TMath::ATan2(y,x);
172 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
173 Double_t r=TMath::Sqrt(x*x + y*y);
175 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
176 r += TMath::Sqrt(x*x + y*y);
177 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
178 r += TMath::Sqrt(x*x + y*y);
179 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
180 r += TMath::Sqrt(x*x + y*y);
183 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
185 for (Int_t j=1; j<nlad+1; j++) {
186 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
187 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
188 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
190 Double_t txyz[3]={0.};
191 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
192 m.LocalToMaster(txyz,xyz);
193 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
194 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
196 if (phi<0) phi+=TMath::TwoPi();
197 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
199 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
200 new(&det) AliITSdetector(r,phi);
201 // compute the real radius (with misalignment)
202 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
204 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
205 mmisal.LocalToMaster(txyz,xyz);
206 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
207 det.SetRmisal(rmisal);
209 } // end loop on detectors
210 } // end loop on ladders
211 fForceSkippingOfLayer[i-1] = 0;
212 } // end loop on layers
215 fI=AliITSgeomTGeo::GetNLayers();
218 fConstraint[0]=1; fConstraint[1]=0;
220 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
221 AliITSReconstructor::GetRecoParam()->GetYVdef(),
222 AliITSReconstructor::GetRecoParam()->GetZVdef()};
223 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
224 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
225 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
226 SetVertex(xyzVtx,ersVtx);
228 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
229 for (Int_t i=0;i<100000;i++){
230 fBestTrackIndex[i]=0;
233 // store positions of centre of SPD modules (in z)
234 // The convetion is that fSPDdetzcentre is ordered from -z to +z
236 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
237 if (tr[2]<0) { // old geom (up to v5asymmPPR)
238 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
239 fSPDdetzcentre[0] = tr[2];
240 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
241 fSPDdetzcentre[1] = tr[2];
242 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
243 fSPDdetzcentre[2] = tr[2];
244 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
245 fSPDdetzcentre[3] = tr[2];
246 } else { // new geom (from v11Hybrid)
247 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
248 fSPDdetzcentre[0] = tr[2];
249 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
250 fSPDdetzcentre[1] = tr[2];
251 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
252 fSPDdetzcentre[2] = tr[2];
253 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
254 fSPDdetzcentre[3] = tr[2];
257 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
258 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
259 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
263 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
264 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
266 if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
267 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
269 // only for plane efficiency evaluation
270 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
271 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
272 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
273 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)!=1)
274 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
276 fPlaneEff = new AliITSPlaneEffSPD();
277 fSPDChipIntPlaneEff = new Bool_t[AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip];
278 for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
280 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
281 else fPlaneEff = new AliITSPlaneEffSSD();
282 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
283 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
284 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
288 fSelectBestMIP03 = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
289 fFlagFakes = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
290 fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
292 fITSPid=new AliITSPIDResponse();
295 //------------------------------------------------------------------------
296 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
298 fBestTrack(tracker.fBestTrack),
299 fTrackToFollow(tracker.fTrackToFollow),
300 fTrackHypothesys(tracker.fTrackHypothesys),
301 fBestHypothesys(tracker.fBestHypothesys),
302 fOriginal(tracker.fOriginal),
303 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
304 fPass(tracker.fPass),
305 fAfterV0(tracker.fAfterV0),
306 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
307 fCoefficients(tracker.fCoefficients),
309 fTrackingPhase(tracker.fTrackingPhase),
310 fUseTGeo(tracker.fUseTGeo),
311 fNtracks(tracker.fNtracks),
312 fFlagFakes(tracker.fFlagFakes),
313 fSelectBestMIP03(tracker.fSelectBestMIP03),
314 fxOverX0Pipe(tracker.fxOverX0Pipe),
315 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
317 fxTimesRhoPipeTrks(0),
318 fxOverX0ShieldTrks(0),
319 fxTimesRhoShieldTrks(0),
320 fxOverX0LayerTrks(0),
321 fxTimesRhoLayerTrks(0),
322 fDebugStreamer(tracker.fDebugStreamer),
323 fITSChannelStatus(tracker.fITSChannelStatus),
324 fkDetTypeRec(tracker.fkDetTypeRec),
325 fPlaneEff(tracker.fPlaneEff) {
327 fOriginal.SetOwner();
330 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
333 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
334 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
337 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
338 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
342 //------------------------------------------------------------------------
343 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
344 //Assignment operator
345 this->~AliITStrackerMI();
346 new(this) AliITStrackerMI(tracker);
350 //------------------------------------------------------------------------
351 AliITStrackerMI::~AliITStrackerMI()
356 if (fCoefficients) delete [] fCoefficients;
357 DeleteTrksMaterialLUT();
358 if (fDebugStreamer) {
359 //fDebugStreamer->Close();
360 delete fDebugStreamer;
362 if(fITSChannelStatus) delete fITSChannelStatus;
363 if(fPlaneEff) delete fPlaneEff;
364 if(fITSPid) delete fITSPid;
365 if (fSPDChipIntPlaneEff) delete [] fSPDChipIntPlaneEff;
368 //------------------------------------------------------------------------
369 void AliITStrackerMI::ReadBadFromDetTypeRec() {
370 //--------------------------------------------------------------------
371 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
373 //--------------------------------------------------------------------
375 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
377 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
379 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
382 if(fITSChannelStatus) delete fITSChannelStatus;
383 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
385 // ITS detectors and chips
386 Int_t i=0,j=0,k=0,ndet=0;
387 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
388 Int_t nBadDetsPerLayer=0;
389 ndet=AliITSgeomTGeo::GetNDetectors(i);
390 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
391 for (k=1; k<ndet+1; k++) {
392 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
393 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
394 if(det.IsBad()) {nBadDetsPerLayer++;}
395 } // end loop on detectors
396 } // end loop on ladders
397 AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
398 } // end loop on layers
402 //------------------------------------------------------------------------
403 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
404 //--------------------------------------------------------------------
405 //This function loads ITS clusters
406 //--------------------------------------------------------------------
408 TClonesArray *clusters = NULL;
409 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
410 clusters=rpcont->FetchClusters(0,cTree);
411 if(!clusters) return 1;
413 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
414 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
417 Int_t i=0,j=0,ndet=0;
419 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
420 ndet=fgLayers[i].GetNdetectors();
421 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
422 for (; j<jmax; j++) {
423 // if (!cTree->GetEvent(j)) continue;
424 clusters = rpcont->UncheckedGetClusters(j);
425 if(!clusters)continue;
426 Int_t ncl=clusters->GetEntriesFast();
427 SignDeltas(clusters,GetZ());
430 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
431 detector=c->GetDetectorIndex();
433 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
435 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
437 AliWarning(Form("Too many clusters on layer %d!",i));
442 // add dead zone "virtual" cluster in SPD, if there is a cluster within
443 // zwindow cm from the dead zone
444 // This method assumes that fSPDdetzcentre is ordered from -z to +z
445 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
446 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
447 Int_t lab[4] = {0,0,0,detector};
448 Int_t info[3] = {0,0,i};
449 Float_t q = 0.; // this identifies virtual clusters
450 Float_t hit[6] = {xdead,
452 static_cast<Float_t>(AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2()),
453 static_cast<Float_t>(AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2()),
456 Bool_t local = kTRUE;
457 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
458 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
459 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
460 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
461 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
462 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
463 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
464 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
465 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
466 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
467 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
468 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
469 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
470 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
471 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
472 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
473 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
474 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
475 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
477 } // "virtual" clusters in SPD
481 fgLayers[i].ResetRoad(); //road defined by the cluster density
482 fgLayers[i].SortClusters();
485 // check whether we have to skip some layers
486 SetForceSkippingOfLayer();
490 //------------------------------------------------------------------------
491 void AliITStrackerMI::UnloadClusters() {
492 //--------------------------------------------------------------------
493 //This function unloads ITS clusters
494 //--------------------------------------------------------------------
495 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
497 //------------------------------------------------------------------------
498 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
499 //--------------------------------------------------------------------
500 // Publishes all pointers to clusters known to the tracker into the
501 // passed object array.
502 // The ownership is not transfered - the caller is not expected to delete
504 //--------------------------------------------------------------------
506 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
507 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
508 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
515 //------------------------------------------------------------------------
516 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
517 //--------------------------------------------------------------------
518 // Correction for the material between the TPC and the ITS
519 //--------------------------------------------------------------------
520 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
521 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
522 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
523 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
524 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
525 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
526 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
527 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
529 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
535 //------------------------------------------------------------------------
536 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
537 //--------------------------------------------------------------------
538 // This functions reconstructs ITS tracks
539 // The clusters must be already loaded !
540 //--------------------------------------------------------------------
542 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
544 fTrackingPhase="Clusters2Tracks";
547 fSelectBestMIP03 = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
548 fFlagFakes = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
549 fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
551 TObjArray itsTracks(15000);
553 fEsd = event; // store pointer to the esd
555 // temporary (for cosmics)
556 if(event->GetVertex()) {
557 TString title = event->GetVertex()->GetTitle();
558 if(title.Contains("cosmics")) {
559 Double_t xyz[3]={GetX(),GetY(),GetZ()};
560 Double_t exyz[3]={0.1,0.1,0.1};
566 {/* Read ESD tracks */
567 Int_t nentr=event->GetNumberOfTracks();
569 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
571 AliESDtrack *esd=event->GetTrack(nentr);
572 // ---- for debugging:
573 //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); }
575 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
576 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
577 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
578 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
579 AliITStrackMI *t = new AliITStrackMI(*esd);
580 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
581 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
584 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
586 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
587 //track - can be V0 according to TPC
589 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
593 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
597 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
601 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
606 t->SetReconstructed(kFALSE);
607 itsTracks.AddLast(t);
608 fOriginal.AddLast(t);
610 } /* End Read ESD tracks */
614 Int_t nentr=itsTracks.GetEntriesFast();
615 fTrackHypothesys.Expand(nentr);
616 fBestHypothesys.Expand(nentr);
617 MakeCoefficients(nentr);
618 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
620 // THE TWO TRACKING PASSES
621 for (fPass=0; fPass<2; fPass++) {
622 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
623 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
624 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
625 if (t==0) continue; //this track has been already tracked
626 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
627 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
628 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
629 if (fConstraint[fPass]) {
630 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
631 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
634 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
635 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
637 ResetTrackToFollow(*t);
640 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
643 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
645 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
646 if (!besttrack) continue;
647 besttrack->SetLabel(tpcLabel);
648 // besttrack->CookdEdx();
650 besttrack->SetFakeRatio(1.);
651 CookLabel(besttrack,0.); //For comparison only
652 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
653 t->SetWinner(besttrack);
655 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
657 t->SetReconstructed(kTRUE);
659 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
661 GetBestHypothesysMIP(itsTracks);
662 } // end loop on the two tracking passes
664 if (fFlagFakes) FlagFakes(itsTracks);
666 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
667 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
672 Int_t entries = fTrackHypothesys.GetEntriesFast();
673 for (Int_t ientry=0; ientry<entries; ientry++) {
674 TObjArray * array =(TObjArray*)fTrackHypothesys.At(ientry);
675 if (array) array->Delete();
676 delete fTrackHypothesys.RemoveAt(ientry);
679 fTrackHypothesys.Delete();
680 entries = fBestHypothesys.GetEntriesFast();
681 for (Int_t ientry=0; ientry<entries; ientry++) {
682 TObjArray * array =(TObjArray*)fBestHypothesys.At(ientry);
683 if (array) array->Delete();
684 delete fBestHypothesys.RemoveAt(ientry);
686 fBestHypothesys.Delete();
688 delete [] fCoefficients;
690 DeleteTrksMaterialLUT();
692 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
694 fTrackingPhase="Default";
698 //------------------------------------------------------------------------
699 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
700 //--------------------------------------------------------------------
701 // This functions propagates reconstructed ITS tracks back
702 // The clusters must be loaded !
703 //--------------------------------------------------------------------
704 fTrackingPhase="PropagateBack";
705 Int_t nentr=event->GetNumberOfTracks();
706 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
707 double bz0 = GetBz();
708 const double kWatchStep=10.; // for larger steps watch arc vs segment difference
711 for (Int_t i=0; i<nentr; i++) {
712 AliESDtrack *esd=event->GetTrack(i);
714 // Start time integral and add distance from current position to vertex
715 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
716 AliITStrackMI t(*esd);
717 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
721 double dxs = xyzTrk[0] - xyzVtx[0];
722 double dys = xyzTrk[1] - xyzVtx[1];
723 double dzs = xyzTrk[2] - xyzVtx[2];
724 // RS: for large segment steps d use approximation of cicrular arc s by
725 // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
726 // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
727 // Hence s^2/d^2 = (1+1/6 p^2)^2
728 dst2 = dxs*dxs + dys*dys;
729 if (dst2 > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
730 double crv = TMath::Abs(esd->GetC(bz0));
731 double fcarc = 1.+crv*crv*dst2/6.;
736 t.StartTimeIntegral();
737 t.AddTimeStep(TMath::Sqrt(dst2));
739 // transfer the time integral to ESD track
740 esd->SetStatus(AliESDtrack::kTIME);
741 Double_t times[AliPID::kSPECIESC];
742 t.GetIntegratedTimes(times,AliPID::kSPECIESC);
743 esd->SetIntegratedTimes(times);
744 esd->SetIntegratedLength(t.GetIntegratedLength());
746 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
748 t.SetExpQ(TMath::Max(0.8*t.GetESDtrack()->GetTPCsignal(),30.));
749 ResetTrackToFollow(t);
751 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
752 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,&t)) {
753 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) continue;
754 // fTrackToFollow.SetLabel(t.GetLabel()); // why do we neet this
755 //fTrackToFollow.CookdEdx();
756 CookLabel(&fTrackToFollow,0.); //For comparison only // why do we need this?
757 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
758 //UseClusters(&fTrackToFollow);
763 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
765 fTrackingPhase="Default";
769 //------------------------------------------------------------------------
770 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
771 //--------------------------------------------------------------------
772 // This functions refits ITS tracks using the
773 // "inward propagated" TPC tracks
774 // The clusters must be loaded !
775 //--------------------------------------------------------------------
776 fTrackingPhase="RefitInward";
778 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
780 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
781 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
783 Int_t nentr=event->GetNumberOfTracks();
784 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
786 // only for PlaneEff and in case of SPD (for FO studies)
787 if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
788 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 &&
789 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) {
790 for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
794 for (Int_t i=0; i<nentr; i++) {
795 AliESDtrack *esd=event->GetTrack(i);
797 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
798 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
799 if (esd->GetStatus()&AliESDtrack::kTPCout)
800 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
802 AliITStrackMI *t = new AliITStrackMI(*esd);
804 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
805 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
810 ResetTrackToFollow(*t);
811 fTrackToFollow.ResetClusters();
813 // ITS standalone tracks
814 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
815 fTrackToFollow.ResetCovariance(10.);
816 // protection for loopers that can have parameters screwed up
817 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
818 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
825 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
826 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
828 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
829 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
830 AliDebug(2," refit OK");
831 fTrackToFollow.SetLabel(t->GetLabel());
832 // fTrackToFollow.CookdEdx();
833 CookdEdx(&fTrackToFollow);
835 CookLabel(&fTrackToFollow,0.0); //For comparison only // RS why do we need this?
838 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
839 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
840 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
841 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
842 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
843 Double_t r[3]={0.,0.,0.};
845 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
852 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
854 fTrackingPhase="Default";
858 //------------------------------------------------------------------------
859 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
860 //--------------------------------------------------------------------
861 // Return pointer to a given cluster
862 //--------------------------------------------------------------------
863 Int_t l=(index & 0xf0000000) >> 28;
864 Int_t c=(index & 0x0fffffff) >> 00;
865 return fgLayers[l].GetCluster(c);
867 //------------------------------------------------------------------------
868 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
869 //--------------------------------------------------------------------
870 // Get track space point with index i
871 //--------------------------------------------------------------------
873 Int_t l=(index & 0xf0000000) >> 28;
874 Int_t c=(index & 0x0fffffff) >> 00;
875 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
876 Int_t idet = cl->GetDetectorIndex();
880 cl->GetGlobalXYZ(xyz);
881 cl->GetGlobalCov(cov);
883 p.SetCharge(cl->GetQ());
884 p.SetDriftTime(cl->GetDriftTime());
885 p.SetChargeRatio(cl->GetChargeRatio());
886 p.SetClusterType(cl->GetClusterType());
887 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
890 iLayer = AliGeomManager::kSPD1;
893 iLayer = AliGeomManager::kSPD2;
896 iLayer = AliGeomManager::kSDD1;
899 iLayer = AliGeomManager::kSDD2;
902 iLayer = AliGeomManager::kSSD1;
905 iLayer = AliGeomManager::kSSD2;
908 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
911 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
912 p.SetVolumeID((UShort_t)volid);
915 //------------------------------------------------------------------------
916 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
917 AliTrackPoint& p, const AliESDtrack *t) {
918 //--------------------------------------------------------------------
919 // Get track space point with index i
920 // (assign error estimated during the tracking)
921 //--------------------------------------------------------------------
923 Int_t l=(index & 0xf0000000) >> 28;
924 Int_t c=(index & 0x0fffffff) >> 00;
925 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
926 Int_t idet = cl->GetDetectorIndex();
928 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
930 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
932 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
933 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
934 Double_t alpha = t->GetAlpha();
935 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
936 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
937 phi += alpha-det.GetPhi();
938 Float_t tgphi = TMath::Tan(phi);
940 Float_t tgl = t->GetTgl(); // tgl about const along track
941 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
943 Float_t errtrky,errtrkz,covyz;
944 Bool_t addMisalErr=kFALSE;
945 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
949 cl->GetGlobalXYZ(xyz);
950 // cl->GetGlobalCov(cov);
951 Float_t pos[3] = {0.,0.,0.};
952 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
953 tmpcl.GetGlobalCov(cov);
956 p.SetCharge(cl->GetQ());
957 p.SetDriftTime(cl->GetDriftTime());
958 p.SetChargeRatio(cl->GetChargeRatio());
959 p.SetClusterType(cl->GetClusterType());
961 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
964 iLayer = AliGeomManager::kSPD1;
967 iLayer = AliGeomManager::kSPD2;
970 iLayer = AliGeomManager::kSDD1;
973 iLayer = AliGeomManager::kSDD2;
976 iLayer = AliGeomManager::kSSD1;
979 iLayer = AliGeomManager::kSSD2;
982 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
985 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
987 p.SetVolumeID((UShort_t)volid);
990 //------------------------------------------------------------------------
991 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
993 //--------------------------------------------------------------------
994 // Follow prolongation tree
995 //--------------------------------------------------------------------
997 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
998 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
1001 AliESDtrack * esd = otrack->GetESDtrack();
1002 if (esd->GetV0Index(0)>0) {
1003 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
1004 // mapping of ESD track is different as ITS track in Containers
1005 // Need something more stable
1006 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
1007 for (Int_t i=0;i<3;i++){
1008 Int_t index = esd->GetV0Index(i);
1009 if (index==0) break;
1010 AliESDv0 * vertex = fEsd->GetV0(index);
1011 if (vertex->GetStatus()<0) continue; // rejected V0
1013 if (esd->GetSign()>0) {
1014 vertex->SetIndex(0,esdindex);
1016 vertex->SetIndex(1,esdindex);
1020 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1022 bestarray = new TObjArray(5);
1023 bestarray->SetOwner();
1024 fBestHypothesys.AddAt(bestarray,esdindex);
1028 //setup tree of the prolongations
1030 const int kMaxTr = 100; //RS
1031 static AliITStrackMI tracks[7][kMaxTr];
1032 AliITStrackMI *currenttrack;
1033 static AliITStrackMI currenttrack1;
1034 static AliITStrackMI currenttrack2;
1035 static AliITStrackMI backuptrack;
1037 Int_t nindexes[7][kMaxTr];
1038 Float_t normalizedchi2[kMaxTr];
1039 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1040 otrack->SetNSkipped(0);
1041 new (&(tracks[6][0])) AliITStrackMI(*otrack);
1043 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1044 Int_t modstatus = 1; // found
1048 // follow prolongations
1049 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1050 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1053 AliITSlayer &layer=fgLayers[ilayer];
1054 Double_t r = layer.GetR();
1060 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1061 // printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1063 if (ntracks[ilayer]>=kMaxTr) break;
1064 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1065 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1066 if (ntracks[ilayer]>15+ilayer){
1067 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1068 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1071 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1073 // material between SSD and SDD, SDD and SPD
1075 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1077 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1081 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1082 Int_t idet=layer.FindDetectorIndex(phi,z);
1084 Double_t trackGlobXYZ1[3];
1085 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1087 // Get the budget to the primary vertex for the current track being prolonged
1088 Double_t budgetToPrimVertex = 0;
1089 double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1092 if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1093 else budgetToPrimVertex = GetEffectiveThickness();
1095 // check if we allow a prolongation without point
1096 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1098 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1099 // propagate to the layer radius
1100 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1101 if(!vtrack->Propagate(xToGo)) continue;
1102 // apply correction for material of the current layer
1103 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1104 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1105 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1106 vtrack->SetClIndex(ilayer,-1);
1107 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1108 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1109 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1111 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1113 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1114 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1120 // track outside layer acceptance in z
1121 if (idet<0) continue;
1123 //propagate to the intersection with the detector plane
1124 const AliITSdetector &det=layer.GetDetector(idet);
1125 new(¤ttrack2) AliITStrackMI(currenttrack1);
1126 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1127 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1128 currenttrack1.SetDetectorIndex(idet);
1129 currenttrack2.SetDetectorIndex(idet);
1130 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1133 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1135 // road in global (rphi,z) [i.e. in tracking ref. system]
1136 Double_t zmin,zmax,ymin,ymax;
1137 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1139 // select clusters in road
1140 layer.SelectClusters(zmin,zmax,ymin,ymax);
1141 //********************
1143 // Define criteria for track-cluster association
1144 Double_t msz = currenttrack1.GetSigmaZ2() +
1145 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1146 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1147 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1148 Double_t msy = currenttrack1.GetSigmaY2() +
1149 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1150 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1151 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1153 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1154 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1156 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1157 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1159 msz = 1./msz; // 1/RoadZ^2
1160 msy = 1./msy; // 1/RoadY^2
1164 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1166 const AliITSRecPoint *cl=0;
1168 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1169 Bool_t deadzoneSPD=kFALSE;
1170 currenttrack = ¤ttrack1;
1172 // check if the road contains a dead zone
1173 Bool_t noClusters = kFALSE;
1174 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1175 if (noClusters) AliDebug(2,"no clusters in road");
1176 Double_t dz=0.5*(zmax-zmin);
1177 Double_t dy=0.5*(ymax-ymin);
1178 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1179 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1180 // create a prolongation without clusters (check also if there are no clusters in the road)
1183 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1184 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1185 updatetrack->SetClIndex(ilayer,-1);
1187 modstatus = 5; // no cls in road
1188 } else if (dead==1) {
1189 modstatus = 7; // holes in z in SPD
1190 } else if (dead==2 || dead==3 || dead==4) {
1191 modstatus = 2; // dead from OCDB
1193 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1194 // apply correction for material of the current layer
1195 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1196 if (constrain) { // apply vertex constrain
1197 updatetrack->SetConstrain(constrain);
1198 Bool_t isPrim = kTRUE;
1199 if (ilayer<4) { // check that it's close to the vertex
1200 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1201 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1202 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1203 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1204 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1206 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1208 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1209 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1212 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1214 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1215 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1217 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1218 updatetrack->SetDeadZoneProbability(ilayer,1.);
1219 } else if (dead==4) { // at least a single dead channel from OCDB
1220 updatetrack->SetDeadZoneProbability(ilayer,0.);
1227 // loop over clusters in the road
1228 while ((cl=layer.GetNextCluster(clidx))!=0) {
1229 if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters
1230 Bool_t changedet =kFALSE;
1231 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1232 Int_t idetc=cl->GetDetectorIndex();
1234 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1235 // take into account misalignment (bring track to real detector plane)
1236 Double_t xTrOrig = currenttrack->GetX();
1237 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1238 // a first cut on track-cluster distance
1239 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1240 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1241 { // cluster not associated to track
1242 AliDebug(2,"not associated");
1243 // MvL: added here as well
1244 // bring track back to ideal detector plane
1245 currenttrack->Propagate(xTrOrig);
1248 // bring track back to ideal detector plane
1249 if (!currenttrack->Propagate(xTrOrig)) continue;
1250 } else { // have to move track to cluster's detector
1251 const AliITSdetector &detc=layer.GetDetector(idetc);
1252 // a first cut on track-cluster distance
1254 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1255 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1256 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1257 continue; // cluster not associated to track
1259 new (&backuptrack) AliITStrackMI(currenttrack2);
1261 currenttrack =¤ttrack2;
1262 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1263 new (currenttrack) AliITStrackMI(backuptrack);
1267 currenttrack->SetDetectorIndex(idetc);
1268 // Get again the budget to the primary vertex
1269 // for the current track being prolonged, if had to change detector
1270 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1273 // calculate track-clusters chi2
1274 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1276 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1277 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1278 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1279 if (ntracks[ilayer]>=kMaxTr) continue;
1280 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1281 updatetrack->SetClIndex(ilayer,-1);
1282 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1284 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1285 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1286 AliDebug(2,"update failed");
1289 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1290 modstatus = 1; // found
1291 } else { // virtual cluster in dead zone
1292 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1293 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1294 modstatus = 7; // holes in z in SPD
1298 Float_t xlocnewdet,zlocnewdet;
1299 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1300 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1303 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1305 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1307 // apply correction for material of the current layer
1308 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1310 if (constrain) { // apply vertex constrain
1311 updatetrack->SetConstrain(constrain);
1312 Bool_t isPrim = kTRUE;
1313 if (ilayer<4) { // check that it's close to the vertex
1314 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1315 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1316 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1317 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1318 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1320 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1322 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1323 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1325 } //apply vertex constrain
1327 } // create new hypothesis
1329 AliDebug(2,"chi2 too large");
1332 } // loop over possible prolongations
1334 // allow one prolongation without clusters
1335 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1336 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1337 // apply correction for material of the current layer
1338 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1339 vtrack->SetClIndex(ilayer,-1);
1340 modstatus = 3; // skipped
1341 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1342 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1344 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1345 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1347 vtrack->IncrementNSkipped();
1352 } // loop over tracks in layer ilayer+1
1354 //loop over track candidates for the current layer
1360 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1361 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1362 if (normalizedchi2[itrack] <
1363 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1367 if (constrain) { // constrain
1368 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1370 } else { // !constrain
1371 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1376 // sort tracks by increasing normalized chi2
1377 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1378 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1379 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1380 // if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1381 if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9);
1382 } // end loop over layers
1386 // Now select tracks to be kept
1388 Int_t max = constrain ? 20 : 5;
1390 // tracks that reach layer 0 (SPD inner)
1391 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1392 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1393 if (track.GetNumberOfClusters()<2) continue;
1394 if (!constrain && track.GetNormChi2(0) >
1395 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1398 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1401 // tracks that reach layer 1 (SPD outer)
1402 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1403 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1404 if (track.GetNumberOfClusters()<4) continue;
1405 if (!constrain && track.GetNormChi2(1) >
1406 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1407 if (constrain) track.IncrementNSkipped();
1409 track.SetD(0,track.GetD(GetX(),GetY()));
1410 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1411 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1412 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1415 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1418 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1420 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1421 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1422 if (track.GetNumberOfClusters()<3) continue;
1423 if (track.GetNormChi2(2) >
1424 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1425 track.SetD(0,track.GetD(GetX(),GetY()));
1426 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1427 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1428 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1430 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1436 // register best track of each layer - important for V0 finder
1438 for (Int_t ilayer=0;ilayer<5;ilayer++){
1439 if (ntracks[ilayer]==0) continue;
1440 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1441 if (track.GetNumberOfClusters()<1) continue;
1442 CookLabel(&track,0);
1443 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1447 // update TPC V0 information
1449 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1450 Float_t fprimvertex[3]={static_cast<Float_t>(GetX()),static_cast<Float_t>(GetY()),static_cast<Float_t>(GetZ())};
1451 for (Int_t i=0;i<3;i++){
1452 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1453 if (index==0) break;
1454 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1455 if (vertex->GetStatus()<0) continue; // rejected V0
1457 if (otrack->GetSign()>0) {
1458 vertex->SetIndex(0,esdindex);
1461 vertex->SetIndex(1,esdindex);
1463 //find nearest layer with track info
1464 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1465 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1466 Int_t nearest = nearestold;
1467 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1468 if (ntracks[nearest]==0){
1473 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1474 if (nearestold<5&&nearest<5){
1475 Bool_t accept = track.GetNormChi2(nearest)<10;
1477 if (track.GetSign()>0) {
1478 vertex->SetParamP(track);
1479 vertex->Update(fprimvertex);
1480 //vertex->SetIndex(0,track.fESDtrack->GetID());
1481 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1483 vertex->SetParamN(track);
1484 vertex->Update(fprimvertex);
1485 //vertex->SetIndex(1,track.fESDtrack->GetID());
1486 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1488 vertex->SetStatus(vertex->GetStatus()+1);
1490 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1497 //------------------------------------------------------------------------
1498 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1500 //--------------------------------------------------------------------
1503 return fgLayers[layer];
1505 //------------------------------------------------------------------------
1506 AliITStrackerMI::AliITSlayer::AliITSlayer():
1536 //--------------------------------------------------------------------
1537 //default AliITSlayer constructor
1538 //--------------------------------------------------------------------
1539 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1540 fClusterWeight[i]=0;
1541 fClusterTracks[0][i]=-1;
1542 fClusterTracks[1][i]=-1;
1543 fClusterTracks[2][i]=-1;
1544 fClusterTracks[3][i]=-1;
1551 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1552 for (Int_t j1=0; j1<6; j1++) {
1553 fClusters5[j1][j]=0;
1554 fClusterIndex5[j1][j]=-1;
1563 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1564 for (Int_t j1=0; j1<11; j1++) {
1565 fClusters10[j1][j]=0;
1566 fClusterIndex10[j1][j]=-1;
1575 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1576 for (Int_t j1=0; j1<21; j1++) {
1577 fClusters20[j1][j]=0;
1578 fClusterIndex20[j1][j]=-1;
1586 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1591 //------------------------------------------------------------------------
1592 AliITStrackerMI::AliITSlayer::
1593 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1622 //--------------------------------------------------------------------
1623 //main AliITSlayer constructor
1624 //--------------------------------------------------------------------
1625 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1626 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1628 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1629 fClusterWeight[i]=0;
1630 fClusterTracks[0][i]=-1;
1631 fClusterTracks[1][i]=-1;
1632 fClusterTracks[2][i]=-1;
1633 fClusterTracks[3][i]=-1;
1641 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1642 for (Int_t j1=0; j1<6; j1++) {
1643 fClusters5[j1][j]=0;
1644 fClusterIndex5[j1][j]=-1;
1653 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1654 for (Int_t j1=0; j1<11; j1++) {
1655 fClusters10[j1][j]=0;
1656 fClusterIndex10[j1][j]=-1;
1665 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1666 for (Int_t j1=0; j1<21; j1++) {
1667 fClusters20[j1][j]=0;
1668 fClusterIndex20[j1][j]=-1;
1676 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1682 //------------------------------------------------------------------------
1683 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1685 fPhiOffset(layer.fPhiOffset),
1686 fNladders(layer.fNladders),
1687 fZOffset(layer.fZOffset),
1688 fNdetectors(layer.fNdetectors),
1689 fDetectors(layer.fDetectors),
1694 fClustersCs(layer.fClustersCs),
1695 fClusterIndexCs(layer.fClusterIndexCs),
1699 fCurrentSlice(layer.fCurrentSlice),
1707 fAccepted(layer.fAccepted),
1709 fMaxSigmaClY(layer.fMaxSigmaClY),
1710 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1711 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1716 //------------------------------------------------------------------------
1717 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1718 //--------------------------------------------------------------------
1719 // AliITSlayer destructor
1720 //--------------------------------------------------------------------
1721 delete [] fDetectors;
1722 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1723 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1724 fClusterWeight[i]=0;
1725 fClusterTracks[0][i]=-1;
1726 fClusterTracks[1][i]=-1;
1727 fClusterTracks[2][i]=-1;
1728 fClusterTracks[3][i]=-1;
1731 //------------------------------------------------------------------------
1732 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1733 //--------------------------------------------------------------------
1734 // This function removes loaded clusters
1735 //--------------------------------------------------------------------
1736 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1737 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1738 fClusterWeight[i]=0;
1739 fClusterTracks[0][i]=-1;
1740 fClusterTracks[1][i]=-1;
1741 fClusterTracks[2][i]=-1;
1742 fClusterTracks[3][i]=-1;
1748 //------------------------------------------------------------------------
1749 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1750 //--------------------------------------------------------------------
1751 // This function reset weights of the clusters
1752 //--------------------------------------------------------------------
1753 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1754 fClusterWeight[i]=0;
1755 fClusterTracks[0][i]=-1;
1756 fClusterTracks[1][i]=-1;
1757 fClusterTracks[2][i]=-1;
1758 fClusterTracks[3][i]=-1;
1760 for (Int_t i=0; i<fN;i++) {
1761 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1762 if (cl&&cl->IsUsed()) cl->Use();
1766 //------------------------------------------------------------------------
1767 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1768 //--------------------------------------------------------------------
1769 // This function calculates the road defined by the cluster density
1770 //--------------------------------------------------------------------
1772 for (Int_t i=0; i<fN; i++) {
1773 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1775 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1777 //------------------------------------------------------------------------
1778 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1779 //--------------------------------------------------------------------
1780 //This function adds a cluster to this layer
1781 //--------------------------------------------------------------------
1782 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1788 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1790 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1791 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1792 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1793 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1794 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1795 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1798 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1799 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1800 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1801 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1805 //------------------------------------------------------------------------
1806 void AliITStrackerMI::AliITSlayer::SortClusters()
1811 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1812 Float_t *z = new Float_t[fN];
1813 Int_t * index = new Int_t[fN];
1815 fMaxSigmaClY=0.; //AD
1816 fMaxSigmaClZ=0.; //AD
1818 for (Int_t i=0;i<fN;i++){
1819 z[i] = fClusters[i]->GetZ();
1820 // save largest errors in y and z for this layer
1821 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1822 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1824 TMath::Sort(fN,z,index,kFALSE);
1825 for (Int_t i=0;i<fN;i++){
1826 clusters[i] = fClusters[index[i]];
1829 for (Int_t i=0;i<fN;i++){
1830 fClusters[i] = clusters[i];
1831 fZ[i] = fClusters[i]->GetZ();
1832 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1833 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1834 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1844 for (Int_t i=0;i<fN;i++){
1845 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1846 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1847 fClusterIndex[i] = i;
1851 fDy5 = (fYB[1]-fYB[0])/5.;
1852 fDy10 = (fYB[1]-fYB[0])/10.;
1853 fDy20 = (fYB[1]-fYB[0])/20.;
1854 for (Int_t i=0;i<6;i++) fN5[i] =0;
1855 for (Int_t i=0;i<11;i++) fN10[i]=0;
1856 for (Int_t i=0;i<21;i++) fN20[i]=0;
1858 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;}
1859 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;}
1860 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;}
1863 for (Int_t i=0;i<fN;i++)
1864 for (Int_t irot=-1;irot<=1;irot++){
1865 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1867 for (Int_t slice=0; slice<6;slice++){
1868 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1869 fClusters5[slice][fN5[slice]] = fClusters[i];
1870 fY5[slice][fN5[slice]] = curY;
1871 fZ5[slice][fN5[slice]] = fZ[i];
1872 fClusterIndex5[slice][fN5[slice]]=i;
1877 for (Int_t slice=0; slice<11;slice++){
1878 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1879 fClusters10[slice][fN10[slice]] = fClusters[i];
1880 fY10[slice][fN10[slice]] = curY;
1881 fZ10[slice][fN10[slice]] = fZ[i];
1882 fClusterIndex10[slice][fN10[slice]]=i;
1887 for (Int_t slice=0; slice<21;slice++){
1888 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1889 fClusters20[slice][fN20[slice]] = fClusters[i];
1890 fY20[slice][fN20[slice]] = curY;
1891 fZ20[slice][fN20[slice]] = fZ[i];
1892 fClusterIndex20[slice][fN20[slice]]=i;
1899 // consistency check
1901 for (Int_t i=0;i<fN-1;i++){
1907 for (Int_t slice=0;slice<21;slice++)
1908 for (Int_t i=0;i<fN20[slice]-1;i++){
1909 if (fZ20[slice][i]>fZ20[slice][i+1]){
1916 //------------------------------------------------------------------------
1917 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1918 //--------------------------------------------------------------------
1919 // This function returns the index of the nearest cluster
1920 //--------------------------------------------------------------------
1923 if (fCurrentSlice<0) {
1932 if (ncl==0) return 0;
1933 Int_t b=0, e=ncl-1, m=(b+e)/2;
1934 for (; b<e; m=(b+e)/2) {
1935 // if (z > fClusters[m]->GetZ()) b=m+1;
1936 if (z > zcl[m]) b=m+1;
1941 //------------------------------------------------------------------------
1942 Bool_t AliITStrackerMI::ComputeRoad(AliITStrackMI* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1943 //--------------------------------------------------------------------
1944 // This function computes the rectangular road for this track
1945 //--------------------------------------------------------------------
1948 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1949 // take into account the misalignment: propagate track to misaligned detector plane
1950 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1952 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1953 TMath::Sqrt(track->GetSigmaZ2() +
1954 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1955 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1956 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1957 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1958 TMath::Sqrt(track->GetSigmaY2() +
1959 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1960 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1961 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1963 // track at boundary between detectors, enlarge road
1964 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1965 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1966 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1967 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1968 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1969 Float_t tgl = TMath::Abs(track->GetTgl());
1970 if (tgl > 1.) tgl=1.;
1971 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1972 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1973 Float_t snp = TMath::Abs(track->GetSnp());
1974 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1975 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1978 // add to the road a term (up to 2-3 mm) to deal with misalignments
1979 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1980 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1982 Double_t r = fgLayers[ilayer].GetR();
1983 zmin = track->GetZ() - dz;
1984 zmax = track->GetZ() + dz;
1985 ymin = track->GetY() + r*det.GetPhi() - dy;
1986 ymax = track->GetY() + r*det.GetPhi() + dy;
1988 // bring track back to idead detector plane
1989 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1993 //------------------------------------------------------------------------
1994 void AliITStrackerMI::AliITSlayer::
1995 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1996 //--------------------------------------------------------------------
1997 // This function sets the "window"
1998 //--------------------------------------------------------------------
2000 Double_t circle=2*TMath::Pi()*fR;
2006 // enlarge road in y by maximum cluster error on this layer (3 sigma)
2007 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
2008 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
2009 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
2010 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
2012 Float_t ymiddle = (fYmin+fYmax)*0.5;
2013 if (ymiddle<fYB[0]) {
2014 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
2015 } else if (ymiddle>fYB[1]) {
2016 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
2022 fClustersCs = fClusters;
2023 fClusterIndexCs = fClusterIndex;
2029 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2030 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2031 if (slice<0) slice=0;
2032 if (slice>20) slice=20;
2033 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2035 fCurrentSlice=slice;
2036 fClustersCs = fClusters20[fCurrentSlice];
2037 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2038 fYcs = fY20[fCurrentSlice];
2039 fZcs = fZ20[fCurrentSlice];
2040 fNcs = fN20[fCurrentSlice];
2045 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2046 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2047 if (slice<0) slice=0;
2048 if (slice>10) slice=10;
2049 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2051 fCurrentSlice=slice;
2052 fClustersCs = fClusters10[fCurrentSlice];
2053 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2054 fYcs = fY10[fCurrentSlice];
2055 fZcs = fZ10[fCurrentSlice];
2056 fNcs = fN10[fCurrentSlice];
2061 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2062 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2063 if (slice<0) slice=0;
2064 if (slice>5) slice=5;
2065 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2067 fCurrentSlice=slice;
2068 fClustersCs = fClusters5[fCurrentSlice];
2069 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2070 fYcs = fY5[fCurrentSlice];
2071 fZcs = fZ5[fCurrentSlice];
2072 fNcs = fN5[fCurrentSlice];
2076 fI = FindClusterIndex(fZmin);
2077 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2083 //------------------------------------------------------------------------
2084 Int_t AliITStrackerMI::AliITSlayer::
2085 FindDetectorIndex(Double_t phi, Double_t z) const {
2086 //--------------------------------------------------------------------
2087 //This function finds the detector crossed by the track
2088 //--------------------------------------------------------------------
2090 if (fZOffset<0) // old geometry
2091 dphi = -(phi-fPhiOffset);
2092 else // new geometry
2093 dphi = phi-fPhiOffset;
2096 if (dphi < 0) dphi += 2*TMath::Pi();
2097 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2098 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2099 if (np>=fNladders) np-=fNladders;
2100 if (np<0) np+=fNladders;
2103 Double_t dz=fZOffset-z;
2104 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2105 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2106 if (nz>=fNdetectors || nz<0) {
2107 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2111 // ad hoc correction for 3rd ladder of SDD inner layer,
2112 // which is reversed (rotated by pi around local y)
2113 // this correction is OK only from AliITSv11Hybrid onwards
2114 if (GetR()>12. && GetR()<20.) { // SDD inner
2115 if(np==2) { // 3rd ladder
2116 Double_t posMod252[3];
2117 AliITSgeomTGeo::GetTranslation(252,posMod252);
2118 // check the Z coordinate of Mod 252: if negative
2119 // (old SDD geometry in AliITSv11Hybrid)
2120 // the swap of numeration whould be applied
2121 if(posMod252[2]<0.){
2122 nz = (fNdetectors-1) - nz;
2126 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2129 return np*fNdetectors + nz;
2131 //------------------------------------------------------------------------
2132 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2134 //--------------------------------------------------------------------
2135 // This function returns clusters within the "window"
2136 //--------------------------------------------------------------------
2138 if (fCurrentSlice<0) {
2139 Double_t rpi2 = 2.*fR*TMath::Pi();
2140 for (Int_t i=fI; i<fImax; i++) {
2143 if (fYmax<y) y -= rpi2;
2144 if (fYmin>y) y += rpi2;
2145 if (y<fYmin) continue;
2146 if (y>fYmax) continue;
2148 // skip clusters that are in "extended" road but they
2149 // 3sigma error does not touch the original road
2150 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2151 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2153 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2156 return fClusters[i];
2159 for (Int_t i=fI; i<fImax; i++) {
2160 if (fYcs[i]<fYmin) continue;
2161 if (fYcs[i]>fYmax) continue;
2162 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2163 ci=fClusterIndexCs[i];
2165 return fClustersCs[i];
2170 //------------------------------------------------------------------------
2171 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2173 //--------------------------------------------------------------------
2174 // This function returns the layer thickness at this point (units X0)
2175 //--------------------------------------------------------------------
2177 x0=AliITSRecoParam::GetX0Air();
2178 if (43<fR&&fR<45) { //SSD2
2181 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2182 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2183 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2184 for (Int_t i=0; i<12; i++) {
2185 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2186 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2190 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2191 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2195 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2196 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2199 if (37<fR&&fR<41) { //SSD1
2202 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2203 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2204 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2205 for (Int_t i=0; i<11; i++) {
2206 if (TMath::Abs(z-3.9*i)<0.15) {
2207 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2211 if (TMath::Abs(z+3.9*i)<0.15) {
2212 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2216 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2217 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2220 if (13<fR&&fR<26) { //SDD
2223 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2225 if (TMath::Abs(y-1.80)<0.55) {
2227 for (Int_t j=0; j<20; j++) {
2228 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2229 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2232 if (TMath::Abs(y+1.80)<0.55) {
2234 for (Int_t j=0; j<20; j++) {
2235 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2236 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2240 for (Int_t i=0; i<4; i++) {
2241 if (TMath::Abs(z-7.3*i)<0.60) {
2243 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2246 if (TMath::Abs(z+7.3*i)<0.60) {
2248 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2253 if (6<fR&&fR<8) { //SPD2
2254 Double_t dd=0.0063; x0=21.5;
2256 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2257 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2259 if (3<fR&&fR<5) { //SPD1
2260 Double_t dd=0.0063; x0=21.5;
2262 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2263 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2268 //------------------------------------------------------------------------
2269 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2271 fRmisal(det.fRmisal),
2273 fSinPhi(det.fSinPhi),
2274 fCosPhi(det.fCosPhi),
2280 fNChips(det.fNChips),
2281 fChipIsBad(det.fChipIsBad)
2285 //------------------------------------------------------------------------
2286 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2287 const AliITSDetTypeRec *detTypeRec)
2289 //--------------------------------------------------------------------
2290 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2291 //--------------------------------------------------------------------
2293 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2294 // while in the tracker they start from 0 for each layer
2295 for(Int_t il=0; il<ilayer; il++)
2296 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2299 if (ilayer==0 || ilayer==1) { // ---------- SPD
2301 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2303 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2306 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2310 // Get calibration from AliITSDetTypeRec
2311 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2312 calib->SetModuleIndex(idet);
2313 AliITSCalibration *calibSPDdead = 0;
2314 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2315 if (calib->IsBad() ||
2316 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2319 // printf("lay %d bad %d\n",ilayer,idet);
2322 // Get segmentation from AliITSDetTypeRec
2323 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2325 // Read info about bad chips
2326 fNChips = segm->GetMaximumChipIndex()+1;
2327 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2328 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2329 fChipIsBad = new Bool_t[fNChips];
2330 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2331 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2332 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2333 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2338 //------------------------------------------------------------------------
2339 Double_t AliITStrackerMI::GetEffectiveThickness()
2341 //--------------------------------------------------------------------
2342 // Returns the thickness between the current layer and the vertex (units X0)
2343 //--------------------------------------------------------------------
2346 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2347 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2348 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2352 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2353 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2357 Double_t xn=fgLayers[fI].GetR();
2358 for (Int_t i=0; i<fI; i++) {
2359 Double_t xi=fgLayers[i].GetR();
2360 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2366 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2367 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2370 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2371 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2376 //------------------------------------------------------------------------
2377 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2379 //--------------------------------------------------------------------
2380 // Returns the array of layers between the current layer and the vertex
2381 //--------------------------------------------------------------------
2384 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2385 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2386 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2391 for (int il=fI;il--;) {
2394 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2395 xMS[nl++] = AliITSRecoParam::GetrInsideShield(1);
2398 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2399 xMS[nl++] = AliITSRecoParam::GetrInsideShield(0);
2402 x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2403 xMS[nl++] = fgLayers[il].GetR();
2408 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2409 xMS[nl++] = AliITSRecoParam::GetrPipe();
2415 //------------------------------------------------------------------------
2416 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2417 //-------------------------------------------------------------------
2418 // This function returns number of clusters within the "window"
2419 //--------------------------------------------------------------------
2421 for (Int_t i=fI; i<fN; i++) {
2422 const AliITSRecPoint *c=fClusters[i];
2423 if (c->GetZ() > fZmax) break;
2424 if (c->IsUsed()) continue;
2425 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2426 Double_t y=fR*det.GetPhi() + c->GetY();
2428 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2429 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2431 if (y<fYmin) continue;
2432 if (y>fYmax) continue;
2437 //------------------------------------------------------------------------
2438 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2439 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2441 //--------------------------------------------------------------------
2442 // This function refits the track "track" at the position "x" using
2443 // the clusters from "clusters"
2444 // If "extra"==kTRUE,
2445 // the clusters from overlapped modules get attached to "track"
2446 // If "planeff"==kTRUE,
2447 // special approach for plane efficiency evaluation is applyed
2448 //--------------------------------------------------------------------
2450 Int_t index[AliITSgeomTGeo::kNLayers];
2452 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2453 Int_t nc=clusters->GetNumberOfClusters();
2454 for (k=0; k<nc; k++) {
2455 Int_t idx=clusters->GetClusterIndex(k);
2456 Int_t ilayer=(idx&0xf0000000)>>28;
2460 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2462 //------------------------------------------------------------------------
2463 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2464 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2466 //--------------------------------------------------------------------
2467 // This function refits the track "track" at the position "x" using
2468 // the clusters from array
2469 // If "extra"==kTRUE,
2470 // the clusters from overlapped modules get attached to "track"
2471 // If "planeff"==kTRUE,
2472 // special approach for plane efficiency evaluation is applyed
2473 //--------------------------------------------------------------------
2474 Int_t index[AliITSgeomTGeo::kNLayers];
2476 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2478 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2479 index[k]=clusters[k];
2482 // special for cosmics and TPC prolonged tracks:
2483 // propagate to the innermost of:
2484 // - innermost layer crossed by the track
2485 // - innermost layer where a cluster was associated to the track
2486 static AliITSRecoParam *repa = NULL;
2488 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2490 repa = AliITSRecoParam::GetHighFluxParam();
2491 AliWarning("Using default AliITSRecoParam class");
2494 Int_t evsp=repa->GetEventSpecie();
2496 if(track->GetESDtrack()) trStatus=track->GetStatus();
2497 Int_t innermostlayer=0;
2498 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2500 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2501 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2502 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2503 index[innermostlayer] >= 0 ) break;
2506 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2509 Int_t modstatus=1; // found
2511 Int_t from, to, step;
2512 if (xx > track->GetX()) {
2513 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2516 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2519 TString dir = (step>0 ? "outward" : "inward");
2521 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2522 AliITSlayer &layer=fgLayers[ilayer];
2523 Double_t r=layer.GetR();
2525 if (step<0 && xx>r) break;
2527 // material between SSD and SDD, SDD and SPD
2528 Double_t hI=ilayer-0.5*step;
2529 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2530 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2531 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2532 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2535 Double_t oldGlobXYZ[3];
2536 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2538 // continue if we are already beyond this layer
2539 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2540 if(step>0 && oldGlobR > r) continue; // going outward
2541 if(step<0 && oldGlobR < r) continue; // going inward
2544 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2546 Int_t idet=layer.FindDetectorIndex(phi,z);
2548 // check if we allow a prolongation without point for large-eta tracks
2549 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2551 modstatus = 4; // out in z
2552 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2553 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2556 // apply correction for material of the current layer
2557 // add time if going outward
2558 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2562 if (idet<0) return kFALSE;
2565 const AliITSdetector &det=layer.GetDetector(idet);
2566 // only for ITS-SA tracks refit
2567 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2569 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2571 track->SetDetectorIndex(idet);
2572 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2574 Double_t dz,zmin,zmax,dy,ymin,ymax;
2576 const AliITSRecPoint *clAcc=0;
2577 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2579 Int_t idx=index[ilayer];
2580 if (idx>=0) { // cluster in this layer
2581 modstatus = 6; // no refit
2582 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2584 if (idet != cl->GetDetectorIndex()) {
2585 idet=cl->GetDetectorIndex();
2586 const AliITSdetector &detc=layer.GetDetector(idet);
2587 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2588 track->SetDetectorIndex(idet);
2589 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2591 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2592 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2596 modstatus = 1; // found
2601 } else { // no cluster in this layer
2603 modstatus = 3; // skipped
2604 // Plane Eff determination:
2605 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2606 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2607 UseTrackForPlaneEff(track,ilayer);
2610 modstatus = 5; // no cls in road
2612 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2613 dz = 0.5*(zmax-zmin);
2614 dy = 0.5*(ymax-ymin);
2615 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2616 if (dead==1) modstatus = 7; // holes in z in SPD
2617 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2622 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2623 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2625 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2628 if (extra && clAcc) { // search for extra clusters in overlapped modules
2629 AliITStrackV2 tmp(*track);
2630 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2631 layer.SelectClusters(zmin,zmax,ymin,ymax);
2633 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2635 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2636 Double_t tolerance=0.1;
2637 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2638 // only clusters in another module! (overlaps)
2639 idetExtra = clExtra->GetDetectorIndex();
2640 if (idet == idetExtra) continue;
2642 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2644 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2645 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2646 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2647 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2649 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2650 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2653 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2654 track->SetExtraModule(ilayer,idetExtra);
2656 } // end search for extra clusters in overlapped modules
2658 // Correct for material of the current layer
2660 // add time if going outward
2661 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2662 track->SetCheckInvariant(kTRUE);
2663 } // end loop on layers
2665 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2669 //------------------------------------------------------------------------
2670 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2673 // calculate normalized chi2
2674 // return NormalizedChi2(track,0);
2677 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2678 // track->fdEdxMismatch=0;
2679 Float_t dedxmismatch =0;
2680 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2682 for (Int_t i = 0;i<6;i++){
2683 if (track->GetClIndex(i)>=0){
2684 Float_t cerry, cerrz;
2685 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2687 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2690 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2691 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2692 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2694 cchi2+=(0.5-ratio)*10.;
2695 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2696 dedxmismatch+=(0.5-ratio)*10.;
2700 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2701 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2702 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2703 if (i<2) chi2+=2*cl->GetDeltaProbability();
2709 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2710 track->SetdEdxMismatch(dedxmismatch);
2714 for (Int_t i = 0;i<4;i++){
2715 if (track->GetClIndex(i)>=0){
2716 Float_t cerry, cerrz;
2717 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2718 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2721 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2722 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2726 for (Int_t i = 4;i<6;i++){
2727 if (track->GetClIndex(i)>=0){
2728 Float_t cerry, cerrz;
2729 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2730 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2733 Float_t cerryb, cerrzb;
2734 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2735 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2738 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2739 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2744 if (track->GetESDtrack()->GetTPCsignal()>85){
2745 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2747 chi2+=(0.5-ratio)*5.;
2750 chi2+=(ratio-2.0)*3;
2754 Double_t match = TMath::Sqrt(track->GetChi22());
2755 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2756 if (!track->GetConstrain()) {
2757 if (track->GetNumberOfClusters()>2) {
2758 match/=track->GetNumberOfClusters()-2.;
2763 if (match<0) match=0;
2765 // penalty factor for missing points (NDeadZone>0), but no penalty
2766 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2767 // or there is a dead from OCDB)
2768 Float_t deadzonefactor = 0.;
2769 if (track->GetNDeadZone()>0.) {
2770 Int_t sumDeadZoneProbability=0;
2771 for(Int_t ilay=0;ilay<6;ilay++) {
2772 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2774 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2775 if(nDeadZoneWithProbNot1>0) {
2776 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2777 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2778 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2780 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2781 deadzonefactor = 3.*(1.1-deadZoneProbability);
2785 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2786 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2787 1./(1.+track->GetNSkipped()));
2788 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2789 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2792 //------------------------------------------------------------------------
2793 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2796 // return matching chi2 between two tracks
2797 Double_t largeChi2=1000.;
2799 AliITStrackMI track3(*track2);
2800 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2802 vec(0,0)=track1->GetY() - track3.GetY();
2803 vec(1,0)=track1->GetZ() - track3.GetZ();
2804 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2805 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2806 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2809 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2810 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2811 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2812 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2813 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2815 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2816 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2817 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2818 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2820 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2821 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2822 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2824 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2825 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2827 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2830 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2831 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2834 //------------------------------------------------------------------------
2835 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2838 // return probability that given point (characterized by z position and error)
2839 // is in SPD dead zone
2840 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2842 Double_t probability = 0.;
2843 Double_t nearestz = 0.,distz=0.;
2844 Int_t nearestzone = -1;
2845 Double_t mindistz = 1000.;
2847 // find closest dead zone
2848 for (Int_t i=0; i<3; i++) {
2849 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2850 if (distz<mindistz) {
2852 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2857 // too far from dead zone
2858 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2861 Double_t zmin, zmax;
2862 if (nearestzone==0) { // dead zone at z = -7
2863 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2864 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2865 } else if (nearestzone==1) { // dead zone at z = 0
2866 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2867 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2868 } else if (nearestzone==2) { // dead zone at z = +7
2869 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2870 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2875 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2877 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2878 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2879 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2882 //------------------------------------------------------------------------
2883 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2886 // calculate normalized chi2
2888 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2890 for (Int_t i = 0;i<6;i++){
2891 if (TMath::Abs(track->GetDy(i))>0){
2892 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2893 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2896 else{chi2[i]=10000;}
2899 TMath::Sort(6,chi2,index,kFALSE);
2900 Float_t max = float(ncl)*fac-1.;
2901 Float_t sumchi=0, sumweight=0;
2902 for (Int_t i=0;i<max+1;i++){
2903 Float_t weight = (i<max)?1.:(max+1.-i);
2904 sumchi+=weight*chi2[index[i]];
2907 Double_t normchi2 = sumchi/sumweight;
2910 //------------------------------------------------------------------------
2911 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2914 // calculate normalized chi2
2915 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2918 for (Int_t i=0;i<6;i++){
2919 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2920 Double_t sy1 = forwardtrack->GetSigmaY(i);
2921 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2922 Double_t sy2 = backtrack->GetSigmaY(i);
2923 Double_t sz2 = backtrack->GetSigmaZ(i);
2924 if (i<2){ sy2=1000.;sz2=1000;}
2926 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2927 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2929 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2930 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2932 res+= nz0*nz0+ny0*ny0;
2935 if (npoints>1) return
2936 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2937 //2*forwardtrack->fNUsed+
2938 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2939 1./(1.+forwardtrack->GetNSkipped()));
2942 //------------------------------------------------------------------------
2943 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2944 //--------------------------------------------------------------------
2945 // Return pointer to a given cluster
2946 //--------------------------------------------------------------------
2947 Int_t l=(index & 0xf0000000) >> 28;
2948 Int_t c=(index & 0x0fffffff) >> 00;
2949 return fgLayers[l].GetWeight(c);
2951 //------------------------------------------------------------------------
2952 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2954 //---------------------------------------------
2955 // register track to the list
2957 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2960 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2961 Int_t index = track->GetClusterIndex(icluster);
2962 Int_t l=(index & 0xf0000000) >> 28;
2963 Int_t c=(index & 0x0fffffff) >> 00;
2964 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2965 for (Int_t itrack=0;itrack<4;itrack++){
2966 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2967 fgLayers[l].SetClusterTracks(itrack,c,id);
2973 //------------------------------------------------------------------------
2974 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2976 //---------------------------------------------
2977 // unregister track from the list
2978 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2979 Int_t index = track->GetClusterIndex(icluster);
2980 Int_t l=(index & 0xf0000000) >> 28;
2981 Int_t c=(index & 0x0fffffff) >> 00;
2982 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2983 for (Int_t itrack=0;itrack<4;itrack++){
2984 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2985 fgLayers[l].SetClusterTracks(itrack,c,-1);
2990 //------------------------------------------------------------------------
2991 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2993 //-------------------------------------------------------------
2994 //get number of shared clusters
2995 //-------------------------------------------------------------
2997 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2998 // mean number of clusters
2999 Float_t *ny = GetNy(id), *nz = GetNz(id);
3002 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
3003 Int_t index = track->GetClusterIndex(icluster);
3004 Int_t l=(index & 0xf0000000) >> 28;
3005 Int_t c=(index & 0x0fffffff) >> 00;
3006 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3007 // if (ny[l]<1.e-13){
3008 // printf("problem\n");
3010 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3014 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3015 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3016 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3018 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3021 deltan = (cl->GetNz()-nz[l]);
3023 if (deltan>2.0) continue; // extended - highly probable shared cluster
3024 weight = 2./TMath::Max(3.+deltan,2.);
3026 for (Int_t itrack=0;itrack<4;itrack++){
3027 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3029 clist[l] = (AliITSRecPoint*)GetCluster(index);
3030 track->SetSharedWeight(l,weight);
3036 track->SetNUsed(shared);
3039 //------------------------------------------------------------------------
3040 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3043 // find first shared track
3045 // mean number of clusters
3046 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
3048 for (Int_t i=0;i<6;i++) overlist[i]=-1;
3049 Int_t sharedtrack=100000;
3050 Int_t tracks[24],trackindex=0;
3051 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3053 for (Int_t icluster=0;icluster<6;icluster++){
3054 if (clusterlist[icluster]<0) continue;
3055 Int_t index = clusterlist[icluster];
3056 Int_t l=(index & 0xf0000000) >> 28;
3057 Int_t c=(index & 0x0fffffff) >> 00;
3058 // if (ny[l]<1.e-13){
3059 // printf("problem\n");
3061 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3062 //if (l>3) continue;
3063 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3066 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3067 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3068 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3070 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3073 deltan = (cl->GetNz()-nz[l]);
3075 if (deltan>2.0) continue; // extended - highly probable shared cluster
3077 for (Int_t itrack=3;itrack>=0;itrack--){
3078 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3079 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3080 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
3085 if (trackindex==0) return -1;
3087 sharedtrack = tracks[0];
3089 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3092 Int_t tracks2[24], cluster[24];
3093 for (Int_t i=0;i<24;i++){ tracks2[i]=-1; cluster[i]=0;}
3096 for (Int_t i=0;i<trackindex;i++){
3097 if (tracks[i]<0) continue;
3098 tracks2[index] = tracks[i];
3100 for (Int_t j=i+1;j<trackindex;j++){
3101 if (tracks[j]<0) continue;
3102 if (tracks[j]==tracks[i]){
3110 for (Int_t i=0;i<index;i++){
3111 if (cluster[index]>max) {
3112 sharedtrack=tracks2[index];
3119 if (sharedtrack>=100000) return -1;
3121 // make list of overlaps
3123 for (Int_t icluster=0;icluster<6;icluster++){
3124 if (clusterlist[icluster]<0) continue;
3125 Int_t index = clusterlist[icluster];
3126 Int_t l=(index & 0xf0000000) >> 28;
3127 Int_t c=(index & 0x0fffffff) >> 00;
3128 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3129 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3131 if (cl->GetNy()>2) continue;
3132 if (cl->GetNz()>2) continue;
3135 if (cl->GetNy()>3) continue;
3136 if (cl->GetNz()>3) continue;
3139 for (Int_t itrack=3;itrack>=0;itrack--){
3140 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3141 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3149 //------------------------------------------------------------------------
3150 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3152 // try to find track hypothesys without conflicts
3153 // with minimal chi2;
3154 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3155 Int_t entries1 = arr1->GetEntriesFast();
3156 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3157 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3158 Int_t entries2 = arr2->GetEntriesFast();
3159 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3161 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3162 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3163 if (track10->Pt()>0.5+track20->Pt()) return track10;
3165 for (Int_t itrack=0;itrack<entries1;itrack++){
3166 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3167 UnRegisterClusterTracks(track,trackID1);
3170 for (Int_t itrack=0;itrack<entries2;itrack++){
3171 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3172 UnRegisterClusterTracks(track,trackID2);
3176 Float_t maxconflicts=6;
3177 Double_t maxchi2 =1000.;
3179 // get weight of hypothesys - using best hypothesy
3182 Int_t list1[6],list2[6];
3183 AliITSRecPoint *clist1[6], *clist2[6] ;
3184 RegisterClusterTracks(track10,trackID1);
3185 RegisterClusterTracks(track20,trackID2);
3186 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3187 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3188 UnRegisterClusterTracks(track10,trackID1);
3189 UnRegisterClusterTracks(track20,trackID2);
3192 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3193 Float_t nerry[6],nerrz[6];
3194 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3195 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3196 for (Int_t i=0;i<6;i++){
3197 if ( (erry1[i]>0) && (erry2[i]>0)) {
3198 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3199 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3201 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3202 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3204 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3205 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3206 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3209 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3210 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3211 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3219 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3220 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3221 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3222 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3224 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3225 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3226 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3228 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3229 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3230 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3233 Double_t sumw = w1+w2;
3237 w1 = (d2+0.5)/(d1+d2+1.);
3238 w2 = (d1+0.5)/(d1+d2+1.);
3240 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3241 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3243 // get pair of "best" hypothesys
3245 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3246 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3248 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3249 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3250 //if (track1->fFakeRatio>0) continue;
3251 RegisterClusterTracks(track1,trackID1);
3252 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3253 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3255 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3256 //if (track2->fFakeRatio>0) continue;
3258 RegisterClusterTracks(track2,trackID2);
3259 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3260 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3261 UnRegisterClusterTracks(track2,trackID2);
3263 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3264 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3265 if (nskipped>0.5) continue;
3267 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3268 if (conflict1+1<cconflict1) continue;
3269 if (conflict2+1<cconflict2) continue;
3273 for (Int_t i=0;i<6;i++){
3275 Float_t c1 =0.; // conflict coeficients
3277 if (clist1[i]&&clist2[i]){
3280 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3283 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3285 c1 = 2./TMath::Max(3.+deltan,2.);
3286 c2 = 2./TMath::Max(3.+deltan,2.);
3292 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3295 deltan = (clist1[i]->GetNz()-nz1[i]);
3297 c1 = 2./TMath::Max(3.+deltan,2.);
3304 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3307 deltan = (clist2[i]->GetNz()-nz2[i]);
3309 c2 = 2./TMath::Max(3.+deltan,2.);
3315 if (TMath::Abs(track1->GetDy(i))>0.) {
3316 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3317 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3318 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3319 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3321 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3324 if (TMath::Abs(track2->GetDy(i))>0.) {
3325 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3326 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3327 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3328 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3331 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3333 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3334 if (chi21>0) sum+=w1;
3335 if (chi22>0) sum+=w2;
3338 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3339 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3340 Double_t normchi2 = 2*conflict+sumchi2/sum;
3341 if ( normchi2 <maxchi2 ){
3344 maxconflicts = conflict;
3348 UnRegisterClusterTracks(track1,trackID1);
3351 // if (maxconflicts<4 && maxchi2<th0){
3352 if (maxchi2<th0*2.){
3353 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3354 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3355 track1->SetChi2MIP(5,maxconflicts);
3356 track1->SetChi2MIP(6,maxchi2);
3357 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3358 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3359 track1->SetChi2MIP(8,index1);
3360 fBestTrackIndex[trackID1] =index1;
3361 UpdateESDtrack(track1, AliESDtrack::kITSin);
3362 original->SetWinner(track1);
3364 else if (track10->GetChi2MIP(0)<th1){
3365 track10->SetChi2MIP(5,maxconflicts);
3366 track10->SetChi2MIP(6,maxchi2);
3367 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3368 UpdateESDtrack(track10,AliESDtrack::kITSin);
3369 original->SetWinner(track10);
3372 for (Int_t itrack=0;itrack<entries1;itrack++){
3373 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3374 UnRegisterClusterTracks(track,trackID1);
3377 for (Int_t itrack=0;itrack<entries2;itrack++){
3378 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3379 UnRegisterClusterTracks(track,trackID2);
3382 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3383 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3384 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3385 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3386 RegisterClusterTracks(track10,trackID1);
3388 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3389 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3390 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3391 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3392 RegisterClusterTracks(track20,trackID2);
3397 //------------------------------------------------------------------------
3398 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3399 //--------------------------------------------------------------------
3400 // This function marks clusters assigned to the track
3401 //--------------------------------------------------------------------
3402 AliTracker::UseClusters(t,from);
3404 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3405 //if (c->GetQ()>2) c->Use();
3406 if (c->GetSigmaZ2()>0.1) c->Use();
3407 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3408 //if (c->GetQ()>2) c->Use();
3409 if (c->GetSigmaZ2()>0.1) c->Use();
3412 //------------------------------------------------------------------------
3413 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3415 //------------------------------------------------------------------
3416 // add track to the list of hypothesys
3417 //------------------------------------------------------------------
3420 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3422 array = new TObjArray(10);
3423 fTrackHypothesys.AddAtAndExpand(array,esdindex);
3425 array->AddLast(track);
3427 //------------------------------------------------------------------------
3428 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3430 //-------------------------------------------------------------------
3431 // compress array of track hypothesys
3432 // keep only maxsize best hypothesys
3433 //-------------------------------------------------------------------
3434 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3435 if (! (fTrackHypothesys.At(esdindex)) ) return;
3436 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3437 Int_t entries = array->GetEntriesFast();
3439 //- find preliminary besttrack as a reference
3440 Float_t minchi2=10000;
3442 AliITStrackMI * besttrack=0;
3444 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3445 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3446 if (!track) continue;
3447 Float_t chi2 = NormalizedChi2(track,0);
3449 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3450 track->SetLabel(tpcLabel);
3452 track->SetFakeRatio(1.);
3453 CookLabel(track,0.); //For comparison only
3455 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3456 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3457 if (track->GetNumberOfClusters()<maxn) continue;
3458 maxn = track->GetNumberOfClusters();
3459 // if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3466 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3467 delete array->RemoveAt(itrack);
3471 if (!besttrack) return;
3474 //take errors of best track as a reference
3475 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3476 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3477 for (Int_t j=0;j<6;j++) {
3478 if (besttrack->GetClIndex(j)>=0){
3479 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3480 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3481 ny[j] = besttrack->GetNy(j);
3482 nz[j] = besttrack->GetNz(j);
3486 // calculate normalized chi2
3488 Float_t * chi2 = new Float_t[entries];
3489 Int_t * index = new Int_t[entries];
3490 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3491 for (Int_t itrack=0;itrack<entries;itrack++){
3492 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3494 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3495 double chi2t = GetNormalizedChi2(track, mode);
3496 track->SetChi2MIP(0,chi2t);
3497 if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3498 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3499 chi2[itrack] = chi2t;
3502 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3503 delete array->RemoveAt(itrack);
3509 TMath::Sort(entries,chi2,index,kFALSE);
3510 besttrack = (AliITStrackMI*)array->At(index[0]);
3511 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3512 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3513 for (Int_t j=0;j<6;j++){
3514 if (besttrack->GetClIndex(j)>=0){
3515 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3516 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3517 ny[j] = besttrack->GetNy(j);
3518 nz[j] = besttrack->GetNz(j);
3523 // calculate one more time with updated normalized errors
3524 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3525 for (Int_t itrack=0;itrack<entries;itrack++){
3526 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3528 double chi2t = GetNormalizedChi2(track, mode);
3529 track->SetChi2MIP(0,chi2t);
3530 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3531 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3532 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3533 chi2[itrack] = chi2t; //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3536 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3537 delete array->RemoveAt(itrack);
3542 entries = array->GetEntriesFast();
3546 TObjArray * newarray = new TObjArray();
3547 TMath::Sort(entries,chi2,index,kFALSE);
3548 besttrack = (AliITStrackMI*)array->At(index[0]);
3550 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3552 for (Int_t j=0;j<6;j++){
3553 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3554 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3555 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3556 ny[j] = besttrack->GetNy(j);
3557 nz[j] = besttrack->GetNz(j);
3560 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3561 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3562 Float_t minn = besttrack->GetNumberOfClusters()-3;
3564 for (Int_t i=0;i<entries;i++){
3565 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3566 if (!track) continue;
3567 if (accepted>maxcut) break;
3568 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3569 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3570 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3571 delete array->RemoveAt(index[i]);
3575 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3576 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3577 if (!shortbest) accepted++;
3579 newarray->AddLast(array->RemoveAt(index[i]));
3580 for (Int_t j=0;j<6;j++){
3582 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3583 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3584 ny[j] = track->GetNy(j);
3585 nz[j] = track->GetNz(j);
3590 delete array->RemoveAt(index[i]);
3594 delete fTrackHypothesys.RemoveAt(esdindex);
3595 fTrackHypothesys.AddAt(newarray,esdindex);
3599 delete fTrackHypothesys.RemoveAt(esdindex);
3605 //------------------------------------------------------------------------
3606 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3608 //-------------------------------------------------------------
3609 // try to find best hypothesy
3610 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3611 // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3612 //-------------------------------------------------------------
3613 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3614 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3615 if (!array) return 0;
3616 Int_t entries = array->GetEntriesFast();
3617 if (!entries) return 0;
3618 Float_t minchi2 = 100000;
3619 AliITStrackMI * besttrack=0;
3621 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3622 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3623 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3624 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3626 for (Int_t i=0;i<entries;i++){
3627 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3628 if (!track) continue;
3629 Float_t sigmarfi,sigmaz;
3630 GetDCASigma(track,sigmarfi,sigmaz);
3631 track->SetDnorm(0,sigmarfi);
3632 track->SetDnorm(1,sigmaz);
3634 track->SetChi2MIP(1,1000000);
3635 track->SetChi2MIP(2,1000000);
3636 track->SetChi2MIP(3,1000000);
3639 backtrack = new(backtrack) AliITStrackMI(*track);
3640 if (track->GetConstrain()) {
3641 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3642 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3643 if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3644 else {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3646 backtrack->ResetCovariance(10.);
3648 backtrack->ResetCovariance(10.);
3650 backtrack->ResetClusters();
3652 Double_t x = original->GetX();
3653 if (!RefitAt(x,backtrack,track)) continue;
3655 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3656 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3657 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3658 track->SetChi22(GetMatchingChi2(backtrack,original));
3660 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3661 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3662 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3665 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3667 //forward track - without constraint
3668 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3669 forwardtrack->ResetClusters();
3671 if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue; // w/o fwd track MIP03 is meaningless
3672 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3673 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3674 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3676 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3677 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3678 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3679 forwardtrack->SetD(0,track->GetD(0));
3680 forwardtrack->SetD(1,track->GetD(1));
3683 AliITSRecPoint* clist[6];
3684 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3685 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3688 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3689 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3690 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3691 track->SetChi2MIP(3,1000);
3694 Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed(); //RS
3695 if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3696 else chi2 += track->GetNUsed();
3698 for (Int_t ichi=0;ichi<5;ichi++){
3699 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3701 if (chi2 < minchi2){
3702 //besttrack = new AliITStrackMI(*forwardtrack);
3704 besttrack->SetLabel(track->GetLabel());
3705 besttrack->SetFakeRatio(track->GetFakeRatio());
3707 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3708 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3709 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3713 delete forwardtrack;
3715 if (!besttrack) return 0;
3718 for (Int_t i=0;i<entries;i++){
3719 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3721 if (!track) continue;
3723 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3724 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3725 // RS: don't apply this cut when fSelectBestMIP03 is on
3726 || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3728 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3729 delete array->RemoveAt(i);
3739 SortTrackHypothesys(esdindex,checkmax,1);
3741 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3742 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3743 besttrack = (AliITStrackMI*)array->At(0);
3744 if (!besttrack) return 0;
3745 besttrack->SetChi2MIP(8,0);
3746 fBestTrackIndex[esdindex]=0;
3747 entries = array->GetEntriesFast();
3748 AliITStrackMI *longtrack =0;
3750 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3751 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3752 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3753 if (!track->GetConstrain()) continue;
3754 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3755 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3756 if (track->GetChi2MIP(0)>4.) continue;
3757 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3760 //if (longtrack) besttrack=longtrack;
3762 // RS do shared cluster analysis here only if the new sharing analysis is not requested
3763 //RRR if (fFlagFakes) return besttrack;
3766 AliITSRecPoint * clist[6];
3767 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3768 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3769 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3770 RegisterClusterTracks(besttrack,esdindex);
3777 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3778 if (sharedtrack>=0){
3780 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3782 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3788 if (shared>2.5) return 0;
3789 if (shared>1.0) return besttrack;
3791 // Don't sign clusters if not gold track
3793 if (!besttrack->IsGoldPrimary()) return besttrack;
3794 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3796 if (fConstraint[fPass]){
3800 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3801 for (Int_t i=0;i<6;i++){
3802 Int_t index = besttrack->GetClIndex(i);
3803 if (index<0) continue;
3804 Int_t ilayer = (index & 0xf0000000) >> 28;
3805 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3806 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3808 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3809 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3810 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3811 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3812 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3813 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3815 Bool_t cansign = kTRUE;
3816 for (Int_t itrack=0;itrack<entries; itrack++){
3817 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3818 if (!track) continue;
3819 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3820 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3826 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3827 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3828 if (!c->IsUsed()) c->Use();
3834 //------------------------------------------------------------------------
3835 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3838 // get "best" hypothesys
3841 Int_t nentries = itsTracks.GetEntriesFast();
3842 for (Int_t i=0;i<nentries;i++){
3843 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3844 if (!track) continue;
3845 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3846 if (!array) continue;
3847 if (array->GetEntriesFast()<=0) continue;
3849 AliITStrackMI* longtrack=0;
3851 Float_t maxchi2=1000;
3852 for (Int_t j=0;j<array->GetEntriesFast();j++){
3853 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3854 if (!trackHyp) continue;
3855 if (trackHyp->GetGoldV0()) {
3856 longtrack = trackHyp; //gold V0 track taken
3859 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3860 Float_t chi2 = trackHyp->GetChi2MIP(0);
3861 if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3862 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3864 if (fAfterV0){ // ??? RS
3865 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3867 if (chi2 > maxchi2) continue;
3868 minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3869 if (fSelectBestMIP03) minn++; // allow next to longest to win
3876 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3877 if (!longtrack) {longtrack = besttrack;}
3878 else besttrack= longtrack;
3882 AliITSRecPoint * clist[6];
3883 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3885 track->SetNUsed(shared);
3886 track->SetNSkipped(besttrack->GetNSkipped());
3887 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3889 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3893 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3894 //if (sharedtrack==-1) sharedtrack=0;
3895 if (sharedtrack>=0) {
3896 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3899 if (besttrack&&fAfterV0) {
3900 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3901 track->SetWinner(besttrack);
3904 if (fConstraint[fPass]) {
3905 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3906 track->SetWinner(besttrack);
3908 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3909 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3910 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3917 //------------------------------------------------------------------------
3918 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3921 // RS: flag those tracks which are suxpected to have fake clusters
3923 const double kThreshPt = 0.5;
3924 AliRefArray *refArr[6];
3926 for (int i=0;i<6;i++) {
3927 int ncl = fgLayers[i].GetNumberOfClusters();
3928 refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3930 Int_t nentries = itsTracks.GetEntriesFast();
3932 // fill cluster->track associations
3933 for (Int_t itr=0;itr<nentries;itr++){
3934 AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3935 if (!track) continue;
3936 AliITStrackMI* trackITS = track->GetWinner();
3937 if (!trackITS) continue;
3938 for (int il=trackITS->GetNumberOfClusters();il--;) {
3939 int idx = trackITS->GetClusterIndex(il);
3940 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3941 // if (c>fgLayers[l].GetNumberOfClusters()) continue;
3942 refArr[l]->AddReference(c, itr);
3946 const UInt_t kMaxRef = 100;
3947 UInt_t crefs[kMaxRef];
3949 // process tracks with shared clusters
3950 for (int itr=0;itr<nentries;itr++){
3951 AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3952 AliITStrackMI* trackH0 = track0->GetWinner();
3953 if (!trackH0) continue;
3954 AliESDtrack* esd0 = track0->GetESDtrack();
3956 for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3957 int idx = trackH0->GetClusterIndex(il);
3958 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3959 ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3960 if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3961 esd0->SetITSSharedFlag(l);
3962 for (int ir=ncrefs;ir--;) {
3963 if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3964 AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3965 AliITStrackMI* trackH1 = track1->GetWinner();
3966 AliESDtrack* esd1 = track1->GetESDtrack();
3967 esd1->SetITSSharedFlag(l);
3969 double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3970 if (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3971 else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3973 // select the one with smallest chi2's product
3974 res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3975 res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3977 if (res<0) esd1->SetITSFakeFlag(); // esd0 is winner
3978 else esd0->SetITSFakeFlag(); // esd1 is winner
3985 for (int i=6;i--;) delete refArr[i];
3990 //------------------------------------------------------------------------
3991 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3992 //--------------------------------------------------------------------
3993 //This function "cooks" a track label. If label<0, this track is fake.
3994 //--------------------------------------------------------------------
3995 const int kMaxLbPerCl = 3;
3996 int lbID[36],lbStat[36];
3997 Int_t nLab=0, nCl = track->GetNumberOfClusters();
3999 // track->SetLabel(-1);
4000 // track->SetFakeRatio(0);
4002 for (Int_t i=0;i<nCl;i++) { // fill all labels
4003 Int_t cindex = track->GetClusterIndex(i);
4004 // Int_t l=(cindex & 0xf0000000) >> 28;
4005 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4007 for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
4008 int trLb = cl->GetLabel(imc);
4010 // search this mc track in already accounted ones
4012 for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
4013 if (iLab<nLab) lbStat[iLab]++;
4018 } // loop over given cluster's labels
4019 } // loop over clusters
4021 if (nLab<1) return; // no labels at all
4024 if (track->GetESDtrack() && track->GetESDtrack()->IsOn(AliESDtrack::kTPCin)){
4025 tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4028 // find majority label
4030 int maxLab=0,tpcLabID=-1;
4031 for (int ilb=nLab;ilb--;) {
4032 int st = lbStat[ilb];
4033 if (lbStat[maxLab]<st) maxLab = ilb;
4034 if (lbID[ilb] == tpcLabel) tpcLabID = ilb;
4036 // if there is an equal choice, prefer ITS label consistent with TPC label
4037 if (tpcLabel>0 && (tpcLabID!=maxLab) && lbStat[maxLab]==lbStat[tpcLabID]) maxLab=tpcLabID;
4039 track->SetFakeRatio(1.-float(lbStat[maxLab])/nCl);
4040 track->SetLabel( lbStat[maxLab]>=nCl-wrong ? lbID[maxLab] : -lbID[maxLab]);
4046 //------------------------------------------------------------------------
4047 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
4048 //--------------------------------------------------------------------
4049 //This function "cooks" a track label. If label<0, this track is fake.
4050 //--------------------------------------------------------------------
4053 if (track->GetESDtrack()){
4054 tpcLabel = track->GetESDtrack()->GetTPCLabel();
4055 ULong_t trStatus=track->GetESDtrack()->GetStatus();
4056 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
4058 track->SetChi2MIP(9,0);
4060 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4061 Int_t cindex = track->GetClusterIndex(i);
4062 Int_t l=(cindex & 0xf0000000) >> 28;
4063 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4065 for (Int_t ind=0;ind<3;ind++){
4066 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
4067 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
4069 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4072 Int_t nclusters = track->GetNumberOfClusters();
4073 if (nclusters > 0) //PH Some tracks don't have any cluster
4074 track->SetFakeRatio(double(nwrong)/double(nclusters));
4075 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4076 track->SetLabel(-tpcLabel);
4078 track->SetLabel(tpcLabel);
4080 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4084 //------------------------------------------------------------------------
4085 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4087 // Fill the dE/dx in this track
4089 track->SetChi2MIP(9,0);
4090 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4091 Int_t cindex = track->GetClusterIndex(i);
4092 Int_t l=(cindex & 0xf0000000) >> 28;
4093 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4094 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4096 for (Int_t ind=0;ind<3;ind++){
4097 if (cl->GetLabel(ind)==lab) isWrong=0;
4099 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4103 track->CookdEdx(low,up);
4105 //------------------------------------------------------------------------
4106 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4108 // Create some arrays
4110 if (fCoefficients) delete []fCoefficients;
4111 fCoefficients = new Float_t[ntracks*48];
4112 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4114 //------------------------------------------------------------------------
4115 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4118 // Compute predicted chi2
4120 // Take into account the mis-alignment (bring track to cluster plane)
4121 Double_t xTrOrig=track->GetX();
4122 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4123 Float_t erry,errz,covyz;
4124 Float_t theta = track->GetTgl();
4125 Float_t phi = track->GetSnp();
4126 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4127 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4128 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()));
4129 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()));
4130 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4131 // Bring the track back to detector plane in ideal geometry
4132 // [mis-alignment will be accounted for in UpdateMI()]
4133 if (!track->Propagate(xTrOrig)) return 1000.;
4135 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4136 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4138 chi2+=0.5*TMath::Min(delta/2,2.);
4139 chi2+=2.*cluster->GetDeltaProbability();
4142 track->SetNy(layer,ny);
4143 track->SetNz(layer,nz);
4144 track->SetSigmaY(layer,erry);
4145 track->SetSigmaZ(layer, errz);
4146 track->SetSigmaYZ(layer,covyz);
4147 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4148 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4152 //------------------------------------------------------------------------
4153 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4158 Int_t layer = (index & 0xf0000000) >> 28;
4159 track->SetClIndex(layer, index);
4160 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4161 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4162 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4163 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4167 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
4170 // Take into account the mis-alignment (bring track to cluster plane)
4171 Double_t xTrOrig=track->GetX();
4172 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4173 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4174 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4175 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
4177 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4180 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4181 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4182 c.SetSigmaYZ(track->GetSigmaYZ(layer));
4185 Int_t updated = track->UpdateMI(&c,chi2,index);
4186 // Bring the track back to detector plane in ideal geometry
4187 if (!track->Propagate(xTrOrig)) return 0;
4189 if(!updated) AliDebug(2,"update failed");
4193 //------------------------------------------------------------------------
4194 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4197 //DCA sigmas parameterization
4198 //to be paramterized using external parameters in future
4201 Double_t curv=track->GetC();
4202 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4203 sigmaz = 0.0110+4.37*TMath::Abs(curv);
4205 //------------------------------------------------------------------------
4206 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4209 // Clusters from delta electrons?
4211 Int_t entries = clusterArray->GetEntriesFast();
4212 if (entries<4) return;
4213 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4214 Int_t layer = cluster->GetLayer();
4215 if (layer>1) return;
4217 Int_t ncandidates=0;
4218 Float_t r = (layer>0)? 7:4;
4220 for (Int_t i=0;i<entries;i++){
4221 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4222 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4223 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4224 index[ncandidates] = i; //candidate to belong to delta electron track
4226 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4227 cl0->SetDeltaProbability(1);
4233 for (Int_t i=0;i<ncandidates;i++){
4234 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4235 if (cl0->GetDeltaProbability()>0.8) continue;
4238 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4239 sumy=sumz=sumy2=sumyz=sumw=0.0;
4240 for (Int_t j=0;j<ncandidates;j++){
4242 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4244 Float_t dz = cl0->GetZ()-cl1->GetZ();
4245 Float_t dy = cl0->GetY()-cl1->GetY();
4246 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4247 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4248 y[ncl] = cl1->GetY();
4249 z[ncl] = cl1->GetZ();
4250 sumy+= y[ncl]*weight;
4251 sumz+= z[ncl]*weight;
4252 sumy2+=y[ncl]*y[ncl]*weight;
4253 sumyz+=y[ncl]*z[ncl]*weight;
4258 if (ncl<4) continue;
4259 Float_t det = sumw*sumy2 - sumy*sumy;
4261 if (TMath::Abs(det)>0.01){
4262 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
4263 Float_t k = (sumyz*sumw - sumy*sumz)/det;
4264 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4267 Float_t z0 = sumyz/sumy;
4268 delta = TMath::Abs(cl0->GetZ()-z0);
4271 cl0->SetDeltaProbability(1-20.*delta);
4275 //------------------------------------------------------------------------
4276 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4281 track->UpdateESDtrack(flags);
4282 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4283 if (oldtrack) delete oldtrack;
4284 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4285 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4286 // printf("Problem\n");
4289 //------------------------------------------------------------------------
4290 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4292 // Get nearest upper layer close to the point xr.
4293 // rough approximation
4295 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4296 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4298 for (Int_t i=0;i<6;i++){
4299 if (radius<kRadiuses[i]){
4306 //------------------------------------------------------------------------
4307 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4308 //--------------------------------------------------------------------
4309 // Fill a look-up table with mean material
4310 //--------------------------------------------------------------------
4314 Double_t point1[3],point2[3];
4315 Double_t phi,cosphi,sinphi,z;
4316 // 0-5 layers, 6 pipe, 7-8 shields
4317 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4318 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4320 Int_t ifirst=0,ilast=0;
4321 if(material.Contains("Pipe")) {
4323 } else if(material.Contains("Shields")) {
4325 } else if(material.Contains("Layers")) {
4328 Error("BuildMaterialLUT","Wrong layer name\n");
4330 const double kAngEps = 1e-4; // tiny slope to avoid tracks strictly normal to Z axis
4331 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4332 Double_t param[5]={0.,0.,0.,0.,0.};
4333 for (Int_t i=0; i<n; i++) {
4334 phi = 2.*TMath::Pi()*gRandom->Rndm();
4335 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4336 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4337 point1[0] = rmin[imat]*cosphi;
4338 point1[1] = rmin[imat]*sinphi;
4340 point2[0] = rmax[imat]*cosphi;
4341 point2[1] = rmax[imat]*sinphi;
4342 point2[2] = z+(rmax[imat]-rmin[imat])*kAngEps;
4343 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4344 if (mparam[1]>999) {n--; continue;} // skip anomalous values in failed propagation
4345 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4347 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4349 fxOverX0Layer[imat] = param[1];
4350 fxTimesRhoLayer[imat] = param[0]*param[4];
4351 } else if(imat==6) {
4352 fxOverX0Pipe = param[1];
4353 fxTimesRhoPipe = param[0]*param[4];
4354 } else if(imat==7) {
4355 fxOverX0Shield[0] = param[1];
4356 fxTimesRhoShield[0] = param[0]*param[4];
4357 } else if(imat==8) {
4358 fxOverX0Shield[1] = param[1];
4359 fxTimesRhoShield[1] = param[0]*param[4];
4363 printf("%s\n",material.Data());
4364 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4365 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4366 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4367 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4368 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4369 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4370 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4371 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4372 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4376 //------------------------------------------------------------------------
4377 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4378 TString direction) {
4379 //-------------------------------------------------------------------
4380 // Propagate beyond beam pipe and correct for material
4381 // (material budget in different ways according to fUseTGeo value)
4382 // Add time if going outward (PropagateTo or PropagateToTGeo)
4383 //-------------------------------------------------------------------
4385 // Define budget mode:
4386 // 0: material from AliITSRecoParam (hard coded)
4387 // 1: material from TGeo in one step (on the fly)
4388 // 2: material from lut
4389 // 3: material from TGeo in one step (same for all hypotheses)
4402 if(fTrackingPhase.Contains("Clusters2Tracks"))
4403 { mode=3; } else { mode=1; }
4406 if(fTrackingPhase.Contains("Clusters2Tracks"))
4407 { mode=3; } else { mode=2; }
4413 if(fTrackingPhase.Contains("Default")) mode=0;
4415 Int_t index=fCurrentEsdTrack;
4417 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4418 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4420 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4422 Double_t xOverX0,x0,lengthTimesMeanDensity;
4426 xOverX0 = AliITSRecoParam::GetdPipe();
4427 x0 = AliITSRecoParam::GetX0Be();
4428 lengthTimesMeanDensity = xOverX0*x0;
4429 lengthTimesMeanDensity *= dir;
4430 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4433 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4436 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4437 xOverX0 = fxOverX0Pipe;
4438 lengthTimesMeanDensity = fxTimesRhoPipe;
4439 lengthTimesMeanDensity *= dir;
4440 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4443 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4444 if(fxOverX0PipeTrks[index]<0) {
4445 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4446 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4447 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4448 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4449 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4452 xOverX0 = fxOverX0PipeTrks[index];
4453 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4454 lengthTimesMeanDensity *= dir;
4455 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4461 //------------------------------------------------------------------------
4462 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4464 TString direction) {
4465 //-------------------------------------------------------------------
4466 // Propagate beyond SPD or SDD shield and correct for material
4467 // (material budget in different ways according to fUseTGeo value)
4468 // Add time if going outward (PropagateTo or PropagateToTGeo)
4469 //-------------------------------------------------------------------
4471 // Define budget mode:
4472 // 0: material from AliITSRecoParam (hard coded)
4473 // 1: material from TGeo in steps of X cm (on the fly)
4474 // X = AliITSRecoParam::GetStepSizeTGeo()
4475 // 2: material from lut
4476 // 3: material from TGeo in one step (same for all hypotheses)
4489 if(fTrackingPhase.Contains("Clusters2Tracks"))
4490 { mode=3; } else { mode=1; }
4493 if(fTrackingPhase.Contains("Clusters2Tracks"))
4494 { mode=3; } else { mode=2; }
4500 if(fTrackingPhase.Contains("Default")) mode=0;
4502 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4504 Int_t shieldindex=0;
4505 if (shield.Contains("SDD")) { // SDDouter
4506 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4508 } else if (shield.Contains("SPD")) { // SPDouter
4509 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4512 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4516 // do nothing if we are already beyond the shield
4517 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4518 if(dir<0 && rTrack > rToGo) return 1; // going outward
4519 if(dir>0 && rTrack < rToGo) return 1; // going inward
4523 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4525 Int_t index=2*fCurrentEsdTrack+shieldindex;
4527 Double_t xOverX0,x0,lengthTimesMeanDensity;
4532 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4533 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4534 lengthTimesMeanDensity = xOverX0*x0;
4535 lengthTimesMeanDensity *= dir;
4536 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4539 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4540 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4543 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4544 xOverX0 = fxOverX0Shield[shieldindex];
4545 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4546 lengthTimesMeanDensity *= dir;
4547 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4550 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4551 if(fxOverX0ShieldTrks[index]<0) {
4552 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4553 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4554 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4555 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4556 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4559 xOverX0 = fxOverX0ShieldTrks[index];
4560 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4561 lengthTimesMeanDensity *= dir;
4562 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4568 //------------------------------------------------------------------------
4569 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4571 Double_t oldGlobXYZ[3],
4572 TString direction) {
4573 //-------------------------------------------------------------------
4574 // Propagate beyond layer and correct for material
4575 // (material budget in different ways according to fUseTGeo value)
4576 // Add time if going outward (PropagateTo or PropagateToTGeo)
4577 //-------------------------------------------------------------------
4579 // Define budget mode:
4580 // 0: material from AliITSRecoParam (hard coded)
4581 // 1: material from TGeo in stepsof X cm (on the fly)
4582 // X = AliITSRecoParam::GetStepSizeTGeo()
4583 // 2: material from lut
4584 // 3: material from TGeo in one step (same for all hypotheses)
4597 if(fTrackingPhase.Contains("Clusters2Tracks"))
4598 { mode=3; } else { mode=1; }
4601 if(fTrackingPhase.Contains("Clusters2Tracks"))
4602 { mode=3; } else { mode=2; }
4608 if(fTrackingPhase.Contains("Default")) mode=0;
4610 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4612 Double_t r=fgLayers[layerindex].GetR();
4613 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4615 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4617 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4619 Int_t index=6*fCurrentEsdTrack+layerindex;
4622 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4625 // back before material (no correction)
4627 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4628 if (!t->GetLocalXat(rOld,xOld)) return 0;
4629 if (!t->Propagate(xOld)) return 0;
4633 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4634 lengthTimesMeanDensity = xOverX0*x0;
4635 lengthTimesMeanDensity *= dir;
4636 // Bring the track beyond the material
4637 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4640 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4641 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4644 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4645 xOverX0 = fxOverX0Layer[layerindex];
4646 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4647 lengthTimesMeanDensity *= dir;
4648 // Bring the track beyond the material
4649 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4652 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4653 if(fxOverX0LayerTrks[index]<0) {
4654 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4655 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4656 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4657 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4658 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4659 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4662 xOverX0 = fxOverX0LayerTrks[index];
4663 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4664 lengthTimesMeanDensity *= dir;
4665 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4672 //------------------------------------------------------------------------
4673 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4674 //-----------------------------------------------------------------
4675 // Initialize LUT for storing material for each prolonged track
4676 //-----------------------------------------------------------------
4677 fxOverX0PipeTrks = new Float_t[ntracks];
4678 fxTimesRhoPipeTrks = new Float_t[ntracks];
4679 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4680 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4681 fxOverX0LayerTrks = new Float_t[ntracks*6];
4682 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4684 for(Int_t i=0; i<ntracks; i++) {
4685 fxOverX0PipeTrks[i] = -1.;
4686 fxTimesRhoPipeTrks[i] = -1.;
4688 for(Int_t j=0; j<ntracks*2; j++) {
4689 fxOverX0ShieldTrks[j] = -1.;
4690 fxTimesRhoShieldTrks[j] = -1.;
4692 for(Int_t k=0; k<ntracks*6; k++) {
4693 fxOverX0LayerTrks[k] = -1.;
4694 fxTimesRhoLayerTrks[k] = -1.;
4701 //------------------------------------------------------------------------
4702 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4703 //-----------------------------------------------------------------
4704 // Delete LUT for storing material for each prolonged track
4705 //-----------------------------------------------------------------
4706 if(fxOverX0PipeTrks) {
4707 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4709 if(fxOverX0ShieldTrks) {
4710 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4713 if(fxOverX0LayerTrks) {
4714 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4716 if(fxTimesRhoPipeTrks) {
4717 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4719 if(fxTimesRhoShieldTrks) {
4720 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4722 if(fxTimesRhoLayerTrks) {
4723 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4727 //------------------------------------------------------------------------
4728 void AliITStrackerMI::SetForceSkippingOfLayer() {
4729 //-----------------------------------------------------------------
4730 // Check if we are forced to skip layers
4731 // either we set to skip them in RecoParam
4732 // or they were off during data-taking
4733 //-----------------------------------------------------------------
4735 const AliEventInfo *eventInfo = GetEventInfo();
4737 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4738 fForceSkippingOfLayer[l] = 0;
4740 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4744 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4745 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4747 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4748 } else if(l==2 || l==3) {
4749 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4751 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4757 //------------------------------------------------------------------------
4758 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4759 Int_t ilayer,Int_t idet) const {
4760 //-----------------------------------------------------------------
4761 // This method is used to decide whether to allow a prolongation
4762 // without clusters, because we want to skip the layer.
4763 // In this case the return value is > 0:
4764 // return 1: the user requested to skip a layer
4765 // return 2: track outside z acceptance
4766 //-----------------------------------------------------------------
4768 if (ForceSkippingOfLayer(ilayer)) return 1;
4770 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4772 if (idet<0 && // out in z
4773 ilayer>innerLayCanSkip &&
4774 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4775 // check if track will cross SPD outer layer
4776 Double_t phiAtSPD2,zAtSPD2;
4777 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4778 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4780 return 2; // always allow skipping, changed on 05.11.2009
4785 //------------------------------------------------------------------------
4786 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4787 Int_t ilayer,Int_t idet,
4788 Double_t dz,Double_t dy,
4789 Bool_t noClusters) const {
4790 //-----------------------------------------------------------------
4791 // This method is used to decide whether to allow a prolongation
4792 // without clusters, because there is a dead zone in the road.
4793 // In this case the return value is > 0:
4794 // return 1: dead zone at z=0,+-7cm in SPD
4795 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4796 // return 2: all road is "bad" (dead or noisy) from the OCDB
4797 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4798 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4799 //-----------------------------------------------------------------
4801 // check dead zones at z=0,+-7cm in the SPD
4802 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4803 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4804 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4805 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4806 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4807 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4808 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4809 for (Int_t i=0; i<3; i++)
4810 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4811 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4812 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4816 // check bad zones from OCDB
4817 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4819 if (idet<0) return 0;
4821 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4824 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4825 if (ilayer==0 || ilayer==1) { // ---------- SPD
4827 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4829 detSizeFactorX *= 2.;
4830 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4833 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4834 if (detType==2) segm->SetLayer(ilayer+1);
4835 Float_t detSizeX = detSizeFactorX*segm->Dx();
4836 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4838 // check if the road overlaps with bad chips
4840 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4841 Float_t zlocmin = zloc-dz;
4842 Float_t zlocmax = zloc+dz;
4843 Float_t xlocmin = xloc-dy;
4844 Float_t xlocmax = xloc+dy;
4845 Int_t chipsInRoad[100];
4847 // check if road goes out of detector
4848 Bool_t touchNeighbourDet=kFALSE;
4849 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4850 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4851 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4852 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4853 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));
4855 // check if this detector is bad
4857 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4858 if(!touchNeighbourDet) {
4859 return 2; // all detectors in road are bad
4861 return 3; // at least one is bad
4865 if(zlocmin>zlocmax)return 0;
4866 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4867 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4868 if (!nChipsInRoad) return 0;
4870 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4871 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4872 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4873 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4874 if (det.IsChipBad(chipsInRoad[iCh])) {
4882 if(!touchNeighbourDet) {
4883 AliDebug(2,"all bad in road");
4884 return 2; // all chips in road are bad
4886 return 3; // at least a bad chip in road
4891 AliDebug(2,"at least a bad in road");
4892 return 3; // at least a bad chip in road
4896 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4897 || !noClusters) return 0;
4899 // There are no clusters in road: check if there is at least
4900 // a bad SPD pixel or SDD anode or SSD strips on both sides
4902 Int_t idetInITS=idet;
4903 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4905 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4906 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4909 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4913 //------------------------------------------------------------------------
4914 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4915 const AliITStrackMI *track,
4916 Float_t &xloc,Float_t &zloc) const {
4917 //-----------------------------------------------------------------
4918 // Gives position of track in local module ref. frame
4919 //-----------------------------------------------------------------
4924 if(idet<0) return kTRUE; // track out of z acceptance of layer
4926 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4928 Int_t lad = Int_t(idet/ndet) + 1;
4930 Int_t det = idet - (lad-1)*ndet + 1;
4932 Double_t xyzGlob[3],xyzLoc[3];
4934 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4935 // take into account the misalignment: xyz at real detector plane
4936 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4938 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4940 xloc = (Float_t)xyzLoc[0];
4941 zloc = (Float_t)xyzLoc[2];
4945 //------------------------------------------------------------------------
4946 //------------------------------------------------------------------------
4947 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer){
4949 // Method to be optimized further:
4950 // Aim: decide whether a track can be used for PlaneEff evaluation
4951 // the decision is taken based on the track quality at the layer under study
4952 // no information on the clusters on this layer has to be used
4953 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4954 // the cut is done on number of sigmas from the boundaries
4956 // Input: Actual track, layer [0,5] under study
4958 // Return: kTRUE if this is a good track
4960 // it will apply a pre-selection to obtain good quality tracks.
4961 // Here also you will have the possibility to put a control on the
4962 // impact point of the track on the basic block, in order to exclude border regions
4963 // this will be done by calling a proper method of the AliITSPlaneEff class.
4965 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4966 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4968 Int_t index[AliITSgeomTGeo::kNLayers];
4970 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4972 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4973 index[k]=clusters[k];
4977 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4978 AliITSlayer &layer=fgLayers[ilayer];
4979 Double_t r=layer.GetR();
4980 AliITStrackMI tmp(*track);
4982 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4983 Int_t nclout=0; Int_t nclin=0;
4984 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4985 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4986 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4987 // if (tmp.GetClIndex(lay)>=0) nclout++;
4988 if(index[lay]>=0)nclout++;
4990 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4991 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4992 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4993 if (index[lay]>=0) nclin++;
4995 Int_t ncl=nclout+nclin;
4996 Bool_t nextout = kFALSE;
4997 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4998 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4999 Bool_t nextin = kFALSE;
5000 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
5001 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
5002 // maximum number of missing clusters allowed in outermost layers
5003 if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
5005 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
5006 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
5008 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
5009 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
5010 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
5011 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
5015 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
5016 Int_t idet=layer.FindDetectorIndex(phi,z);
5017 if(idet<0) { AliInfo(Form("cannot find detector"));
5020 // here check if it has good Chi Square.
5022 //propagate to the intersection with the detector plane
5023 const AliITSdetector &det=layer.GetDetector(idet);
5024 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
5028 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
5029 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5030 if(key>fPlaneEff->Nblock()) return kFALSE;
5031 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
5032 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
5034 // DEFINITION OF SEARCH ROAD FOR accepting a track
5036 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
5037 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
5038 Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXFromBoundaryPlaneEff();
5039 Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZFromBoundaryPlaneEff();
5040 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
5041 // those are precisions in the tracking reference system
5042 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
5043 // Use it also for the module reference system, as it is done for RecPoints
5045 if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigFrmBndPlaneEff()){
5046 if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) dx -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
5049 if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) dz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
5053 // exclude tracks at boundary between detectors
5054 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
5055 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
5056 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
5057 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
5058 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5059 if ( (locx-dx < blockXmn+boundaryWidth) ||
5060 (locx+dx > blockXmx-boundaryWidth) ||
5061 (locz-dz < blockZmn+boundaryWidth) ||
5062 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5065 const AliESDEvent *myesd = ((AliESDtrack*)tmp.GetESDtrack())->GetESDEvent();
5067 if (CorrectForPipeMaterial(&tmp,"inward")) {
5068 const AliESDVertex* vtx = myesd->GetVertex();
5069 if(!vtx) return kFALSE;
5070 Double_t ddz[2],cov[3];
5072 if(!tmp.PropagateToDCA(vtx,AliTracker::GetBz(),maxD,ddz,cov)) return kFALSE;
5073 if(TMath::Abs(ddz[0])>=AliITSReconstructor::GetRecoParam()->GetDCACutPlaneEff()) return kFALSE;
5075 Double_t covar[6]; vtx->GetCovMatrix(covar);
5076 Double_t p[2]={tmp.GetParameter()[0]-ddz[0],tmp.GetParameter()[1]-ddz[1]};
5077 Double_t c[3]={covar[2],0.,covar[5]};
5078 Double_t chi2= ((AliESDtrack*)tmp.GetESDtrack())->GetPredictedChi2(p,c);
5079 if (chi2>AliITSReconstructor::GetRecoParam()->GetVertexChi2CutPlaneEff()) return kFALSE; // Use this to cut on chi^2
5081 } else return kFALSE;
5087 //------------------------------------------------------------------------
5088 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
5090 // This Method has to be optimized! For the time-being it uses the same criteria
5091 // as those used in the search of extra clusters for overlapping modules.
5093 // Method Purpose: estabilish whether a track has produced a recpoint or not
5094 // in the layer under study (For Plane efficiency)
5096 // inputs: AliITStrackMI* track (pointer to a usable track)
5098 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5099 // traversed by this very track. In details:
5100 // - if a cluster can be associated to the track then call
5101 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5102 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5105 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5106 AliITSlayer &layer=fgLayers[ilayer];
5107 Double_t r=layer.GetR();
5108 AliITStrackMI tmp(*track);
5112 if (!tmp.GetPhiZat(r,phi,z)) return;
5113 Int_t idet=layer.FindDetectorIndex(phi,z);
5115 if(idet<0) { AliInfo(Form("cannot find detector"));
5119 //propagate to the intersection with the detector plane
5120 const AliITSdetector &det=layer.GetDetector(idet);
5121 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5125 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5127 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5128 TMath::Sqrt(tmp.GetSigmaZ2() +
5129 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5130 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5131 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5132 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5133 TMath::Sqrt(tmp.GetSigmaY2() +
5134 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5135 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5136 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5138 // road in global (rphi,z) [i.e. in tracking ref. system]
5139 Double_t zmin = tmp.GetZ() - dz;
5140 Double_t zmax = tmp.GetZ() + dz;
5141 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5142 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5144 // select clusters in road
5145 layer.SelectClusters(zmin,zmax,ymin,ymax);
5147 // Define criteria for track-cluster association
5148 Double_t msz = tmp.GetSigmaZ2() +
5149 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5150 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5151 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5152 Double_t msy = tmp.GetSigmaY2() +
5153 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5154 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5155 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5156 if (tmp.GetConstrain()) {
5157 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5158 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5160 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5161 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5164 if(AliITSReconstructor::GetRecoParam()->GetSwitchOffStdSearchClusPlaneEff()){
5165 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXSearchClusterPlaneEff();
5166 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZSearchClusterPlaneEff();
5167 Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXSearchClusterPlaneEff();
5168 Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZSearchClusterPlaneEff();
5169 msy = nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
5170 msz = nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
5172 if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigSrhClusPlaneEff()){
5173 if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) msy -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
5176 if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) msz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
5185 if(msz==0 || msy==0){AliWarning("UseTrackForPlaneEff: null search frame"); return;}
5187 msz = 1./msz; // 1/RoadZ^2
5188 msy = 1./msy; // 1/RoadY^2
5190 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5192 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5193 //Double_t tolerance=0.2;
5194 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5195 idetc = cl->GetDetectorIndex();
5196 if(idet!=idetc) continue;
5197 //Int_t ilay = cl->GetLayer();
5199 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5200 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5202 Double_t chi2=tmp.GetPredictedChi2(cl);
5203 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5207 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5209 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5210 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5211 if(key>fPlaneEff->Nblock()) return;
5212 Bool_t found=kFALSE;
5215 while ((cl=layer.GetNextCluster(clidx))!=0) {
5216 idetc = cl->GetDetectorIndex();
5217 if(idet!=idetc) continue;
5218 // here real control to see whether the cluster can be associated to the track.
5219 // cluster not associated to track
5220 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5221 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5222 // calculate track-clusters chi2
5223 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5224 // in particular, the error associated to the cluster
5225 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5227 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5229 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5230 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5231 // track->SetExtraModule(ilayer,idetExtra);
5233 if(!fPlaneEff->UpDatePlaneEff(found,key))
5234 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5236 // this for FO efficiency studies (only for SPD) //
5237 UInt_t keyFO=999999;
5238 Bool_t foundFO=kFALSE;
5239 if(ilayer<2){ //ONLY SPD layers for FastOr studies
5240 TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5241 Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5242 if(!fSPDChipIntPlaneEff[key]){
5243 AliITSPlaneEffSPD spd;
5244 keyFO = spd.SwitchChipKeyNumbering(key);
5245 if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5246 keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5247 if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5248 AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5251 if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5252 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5258 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5259 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5260 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5261 Int_t cltype[2]={-999,-999};
5264 Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5268 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5269 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5272 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5273 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5274 cltype[0]=layer.GetCluster(ci)->GetNy();
5275 cltype[1]=layer.GetCluster(ci)->GetNz();
5277 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5278 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5279 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5280 // It is computed properly by calling the method
5281 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5283 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5284 //if (tmp.PropagateTo(x,0.,0.)) {
5285 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5286 AliCluster c(*layer.GetCluster(ci));
5287 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5288 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5289 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5290 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5291 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5294 // Compute the angles between the track and the module
5295 // compute the angle "in phi direction", i.e. the angle in the transverse plane
5296 // between the normal to the module and the projection (in the transverse plane) of the
5298 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5299 Float_t tgl = tmp.GetTgl();
5300 Float_t phitr = tmp.GetSnp();
5301 phitr = TMath::ASin(phitr);
5302 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5304 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5305 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5307 alpha = tmp.GetAlpha();
5308 Double_t phiglob = alpha+phitr;
5310 p[0] = TMath::Cos(phiglob);
5311 p[1] = TMath::Sin(phiglob);
5313 TVector3 pvec(p[0],p[1],p[2]);
5314 TVector3 normvec(rot[1],rot[4],rot[7]);
5315 Double_t angle = pvec.Angle(normvec);
5317 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5318 angle *= 180./TMath::Pi();
5321 TVector3 pt(p[0],p[1],0);
5322 TVector3 normt(rot[1],rot[4],0);
5323 Double_t anglet = pt.Angle(normt);
5325 Double_t phiPt = TMath::ATan2(p[1],p[0]);
5326 if(phiPt<0)phiPt+=2.*TMath::Pi();
5327 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5328 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5329 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5330 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
5331 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5332 anglet *= 180./TMath::Pi();
5334 angleModTrack[2]=(Float_t) angle;
5335 angleModTrack[0]=(Float_t) anglet;
5336 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5337 angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5338 angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5339 angleModTrack[1]*=180./TMath::Pi(); // in degree
5341 fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
5343 // For FO efficiency studies of SPD
5344 if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
5346 if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5350 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5352 // find the MC cluster for the label
5353 return fgLayers[lr].FindClusterForLabel(label,store);
5357 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5359 // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5360 strncpy(patt,"......",6);
5362 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5365 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5366 Int_t cindex = track->GetClusterIndex(i);
5367 Int_t l=(cindex & 0xf0000000) >> 28;
5368 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5370 for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5371 patt[l] = isWrong ? 'f':'c';
5377 //------------------------------------------------------------------------
5378 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5380 //--------------------------------------------------------------------
5383 for (int ic=0;ic<fN;ic++) {
5384 const AliITSRecPoint *cl = GetCluster(ic);
5385 for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5387 if (store) store[nfound] = ic;