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"
63 ClassImp(AliITStrackerMI)
65 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
67 AliITStrackerMI::AliITStrackerMI():AliTracker(),
77 fLastLayerToTrackTo(0),
80 fTrackingPhase("Default"),
84 fSelectBestMIP03(kFALSE),
85 fUseImproveKalman(kFALSE),
89 fxTimesRhoPipeTrks(0),
90 fxOverX0ShieldTrks(0),
91 fxTimesRhoShieldTrks(0),
93 fxTimesRhoLayerTrks(0),
98 fSPDChipIntPlaneEff(0),
102 //Default constructor
104 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
106 fxOverX0Shield[i]=-1.;
107 fxTimesRhoShield[i]=-1.;
110 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
111 fOriginal.SetOwner();
112 for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
113 for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
114 fITSPid=new AliITSPIDResponse();
117 //------------------------------------------------------------------------
118 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
119 fI(AliITSgeomTGeo::GetNLayers()),
128 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
131 fTrackingPhase("Default"),
135 fSelectBestMIP03(kFALSE),
136 fUseImproveKalman(kFALSE),
140 fxTimesRhoPipeTrks(0),
141 fxOverX0ShieldTrks(0),
142 fxTimesRhoShieldTrks(0),
143 fxOverX0LayerTrks(0),
144 fxTimesRhoLayerTrks(0),
146 fITSChannelStatus(0),
149 fSPDChipIntPlaneEff(0),
151 //--------------------------------------------------------------------
152 //This is the AliITStrackerMI constructor
153 //--------------------------------------------------------------------
155 AliWarning("\"geom\" is actually a dummy argument !");
158 fOriginal.SetOwner();
162 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
163 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
164 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
166 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
167 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
168 Double_t poff=TMath::ATan2(y,x);
171 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
172 Double_t r=TMath::Sqrt(x*x + y*y);
174 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
175 r += TMath::Sqrt(x*x + y*y);
176 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
177 r += TMath::Sqrt(x*x + y*y);
178 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
179 r += TMath::Sqrt(x*x + y*y);
182 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
184 for (Int_t j=1; j<nlad+1; j++) {
185 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
186 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
187 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
189 Double_t txyz[3]={0.};
190 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
191 m.LocalToMaster(txyz,xyz);
192 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
193 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
195 if (phi<0) phi+=TMath::TwoPi();
196 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
198 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
199 new(&det) AliITSdetector(r,phi);
200 // compute the real radius (with misalignment)
201 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
203 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
204 mmisal.LocalToMaster(txyz,xyz);
205 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
206 det.SetRmisal(rmisal);
208 } // end loop on detectors
209 } // end loop on ladders
210 fForceSkippingOfLayer[i-1] = 0;
211 } // end loop on layers
214 fI=AliITSgeomTGeo::GetNLayers();
217 fConstraint[0]=1; fConstraint[1]=0;
219 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
220 AliITSReconstructor::GetRecoParam()->GetYVdef(),
221 AliITSReconstructor::GetRecoParam()->GetZVdef()};
222 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
223 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
224 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
225 SetVertex(xyzVtx,ersVtx);
227 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
228 for (Int_t i=0;i<100000;i++){
229 fBestTrackIndex[i]=0;
232 // store positions of centre of SPD modules (in z)
233 // The convetion is that fSPDdetzcentre is ordered from -z to +z
235 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
236 if (tr[2]<0) { // old geom (up to v5asymmPPR)
237 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
238 fSPDdetzcentre[0] = tr[2];
239 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
240 fSPDdetzcentre[1] = tr[2];
241 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
242 fSPDdetzcentre[2] = tr[2];
243 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
244 fSPDdetzcentre[3] = tr[2];
245 } else { // new geom (from v11Hybrid)
246 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
247 fSPDdetzcentre[0] = tr[2];
248 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
249 fSPDdetzcentre[1] = tr[2];
250 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
251 fSPDdetzcentre[2] = tr[2];
252 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
253 fSPDdetzcentre[3] = tr[2];
256 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
257 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
258 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
262 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
263 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
265 if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
266 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
268 // only for plane efficiency evaluation
269 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
270 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
271 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
272 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
273 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
275 fPlaneEff = new AliITSPlaneEffSPD();
276 fSPDChipIntPlaneEff = new Bool_t[AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip];
277 for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
279 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
280 else fPlaneEff = new AliITSPlaneEffSSD();
281 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
282 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
283 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
287 fSelectBestMIP03 = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
288 fFlagFakes = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
289 fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
291 fITSPid=new AliITSPIDResponse();
294 //------------------------------------------------------------------------
295 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
297 fBestTrack(tracker.fBestTrack),
298 fTrackToFollow(tracker.fTrackToFollow),
299 fTrackHypothesys(tracker.fTrackHypothesys),
300 fBestHypothesys(tracker.fBestHypothesys),
301 fOriginal(tracker.fOriginal),
302 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
303 fPass(tracker.fPass),
304 fAfterV0(tracker.fAfterV0),
305 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
306 fCoefficients(tracker.fCoefficients),
308 fTrackingPhase(tracker.fTrackingPhase),
309 fUseTGeo(tracker.fUseTGeo),
310 fNtracks(tracker.fNtracks),
311 fFlagFakes(tracker.fFlagFakes),
312 fSelectBestMIP03(tracker.fSelectBestMIP03),
313 fxOverX0Pipe(tracker.fxOverX0Pipe),
314 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
316 fxTimesRhoPipeTrks(0),
317 fxOverX0ShieldTrks(0),
318 fxTimesRhoShieldTrks(0),
319 fxOverX0LayerTrks(0),
320 fxTimesRhoLayerTrks(0),
321 fDebugStreamer(tracker.fDebugStreamer),
322 fITSChannelStatus(tracker.fITSChannelStatus),
323 fkDetTypeRec(tracker.fkDetTypeRec),
324 fPlaneEff(tracker.fPlaneEff) {
326 fOriginal.SetOwner();
329 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
332 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
333 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
336 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
337 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
341 //------------------------------------------------------------------------
342 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
343 //Assignment operator
344 this->~AliITStrackerMI();
345 new(this) AliITStrackerMI(tracker);
349 //------------------------------------------------------------------------
350 AliITStrackerMI::~AliITStrackerMI()
355 if (fCoefficients) delete [] fCoefficients;
356 DeleteTrksMaterialLUT();
357 if (fDebugStreamer) {
358 //fDebugStreamer->Close();
359 delete fDebugStreamer;
361 if(fITSChannelStatus) delete fITSChannelStatus;
362 if(fPlaneEff) delete fPlaneEff;
363 if(fITSPid) delete fITSPid;
364 if (fSPDChipIntPlaneEff) delete [] fSPDChipIntPlaneEff;
367 //------------------------------------------------------------------------
368 void AliITStrackerMI::ReadBadFromDetTypeRec() {
369 //--------------------------------------------------------------------
370 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
372 //--------------------------------------------------------------------
374 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
376 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
378 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
381 if(fITSChannelStatus) delete fITSChannelStatus;
382 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
384 // ITS detectors and chips
385 Int_t i=0,j=0,k=0,ndet=0;
386 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
387 Int_t nBadDetsPerLayer=0;
388 ndet=AliITSgeomTGeo::GetNDetectors(i);
389 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
390 for (k=1; k<ndet+1; k++) {
391 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
392 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
393 if(det.IsBad()) {nBadDetsPerLayer++;}
394 } // end loop on detectors
395 } // end loop on ladders
396 AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
397 } // end loop on layers
401 //------------------------------------------------------------------------
402 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
403 //--------------------------------------------------------------------
404 //This function loads ITS clusters
405 //--------------------------------------------------------------------
407 TClonesArray *clusters = NULL;
408 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
409 clusters=rpcont->FetchClusters(0,cTree);
410 if(!clusters) return 1;
412 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
413 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
416 Int_t i=0,j=0,ndet=0;
418 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
419 ndet=fgLayers[i].GetNdetectors();
420 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
421 for (; j<jmax; j++) {
422 // if (!cTree->GetEvent(j)) continue;
423 clusters = rpcont->UncheckedGetClusters(j);
424 if(!clusters)continue;
425 Int_t ncl=clusters->GetEntriesFast();
426 SignDeltas(clusters,GetZ());
429 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
430 detector=c->GetDetectorIndex();
432 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
434 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
436 AliWarning(Form("Too many clusters on layer %d!",i));
441 // add dead zone "virtual" cluster in SPD, if there is a cluster within
442 // zwindow cm from the dead zone
443 // This method assumes that fSPDdetzcentre is ordered from -z to +z
444 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
445 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
446 Int_t lab[4] = {0,0,0,detector};
447 Int_t info[3] = {0,0,i};
448 Float_t q = 0.; // this identifies virtual clusters
449 Float_t hit[6] = {xdead,
451 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
452 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
455 Bool_t local = kTRUE;
456 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
457 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
458 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
459 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
460 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
461 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
462 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
463 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
464 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
465 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
466 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
467 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
468 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
469 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
470 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
471 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
472 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
473 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
474 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
476 } // "virtual" clusters in SPD
480 fgLayers[i].ResetRoad(); //road defined by the cluster density
481 fgLayers[i].SortClusters();
484 // check whether we have to skip some layers
485 SetForceSkippingOfLayer();
489 //------------------------------------------------------------------------
490 void AliITStrackerMI::UnloadClusters() {
491 //--------------------------------------------------------------------
492 //This function unloads ITS clusters
493 //--------------------------------------------------------------------
494 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
496 //------------------------------------------------------------------------
497 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
498 //--------------------------------------------------------------------
499 // Publishes all pointers to clusters known to the tracker into the
500 // passed object array.
501 // The ownership is not transfered - the caller is not expected to delete
503 //--------------------------------------------------------------------
505 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
506 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
507 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
514 //------------------------------------------------------------------------
515 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
516 //--------------------------------------------------------------------
517 // Correction for the material between the TPC and the ITS
518 //--------------------------------------------------------------------
519 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
520 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
521 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
522 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
523 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
524 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
525 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
526 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
528 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
534 //------------------------------------------------------------------------
535 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
536 //--------------------------------------------------------------------
537 // This functions reconstructs ITS tracks
538 // The clusters must be already loaded !
539 //--------------------------------------------------------------------
541 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
543 fTrackingPhase="Clusters2Tracks";
546 fSelectBestMIP03 = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
547 fFlagFakes = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
548 fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
550 TObjArray itsTracks(15000);
552 fEsd = event; // store pointer to the esd
554 // temporary (for cosmics)
555 if(event->GetVertex()) {
556 TString title = event->GetVertex()->GetTitle();
557 if(title.Contains("cosmics")) {
558 Double_t xyz[3]={GetX(),GetY(),GetZ()};
559 Double_t exyz[3]={0.1,0.1,0.1};
565 {/* Read ESD tracks */
566 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
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 // look at the ESD mass hypothesys !
585 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
587 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
589 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
590 //track - can be V0 according to TPC
592 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
596 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
600 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
604 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
609 t->SetReconstructed(kFALSE);
610 itsTracks.AddLast(t);
611 fOriginal.AddLast(t);
613 } /* End Read ESD tracks */
617 Int_t nentr=itsTracks.GetEntriesFast();
618 fTrackHypothesys.Expand(nentr);
619 fBestHypothesys.Expand(nentr);
620 MakeCoefficients(nentr);
621 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
623 // THE TWO TRACKING PASSES
624 for (fPass=0; fPass<2; fPass++) {
625 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
626 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
627 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
628 if (t==0) continue; //this track has been already tracked
629 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
630 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
631 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
632 if (fConstraint[fPass]) {
633 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
634 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
637 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
638 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
640 ResetTrackToFollow(*t);
643 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
646 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
648 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
649 if (!besttrack) continue;
650 besttrack->SetLabel(tpcLabel);
651 // besttrack->CookdEdx();
653 besttrack->SetFakeRatio(1.);
654 CookLabel(besttrack,0.); //For comparison only
655 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
656 t->SetWinner(besttrack);
658 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
660 t->SetReconstructed(kTRUE);
662 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
664 GetBestHypothesysMIP(itsTracks);
665 } // end loop on the two tracking passes
667 if (fFlagFakes) FlagFakes(itsTracks);
669 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
670 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
675 Int_t entries = fTrackHypothesys.GetEntriesFast();
676 for (Int_t ientry=0; ientry<entries; ientry++) {
677 TObjArray * array =(TObjArray*)fTrackHypothesys.At(ientry);
678 if (array) array->Delete();
679 delete fTrackHypothesys.RemoveAt(ientry);
682 fTrackHypothesys.Delete();
683 entries = fBestHypothesys.GetEntriesFast();
684 for (Int_t ientry=0; ientry<entries; ientry++) {
685 TObjArray * array =(TObjArray*)fBestHypothesys.At(ientry);
686 if (array) array->Delete();
687 delete fBestHypothesys.RemoveAt(ientry);
689 fBestHypothesys.Delete();
691 delete [] fCoefficients;
693 DeleteTrksMaterialLUT();
695 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
697 fTrackingPhase="Default";
701 //------------------------------------------------------------------------
702 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
703 //--------------------------------------------------------------------
704 // This functions propagates reconstructed ITS tracks back
705 // The clusters must be loaded !
706 //--------------------------------------------------------------------
707 fTrackingPhase="PropagateBack";
708 Int_t nentr=event->GetNumberOfTracks();
709 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
712 for (Int_t i=0; i<nentr; i++) {
713 AliESDtrack *esd=event->GetTrack(i);
715 // Start time integral and add distance from current position to vertex
716 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
717 AliITStrackMI *t = new AliITStrackMI(*esd);
718 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
721 for (Int_t icoord=0; icoord<3; icoord++) {Double_t di = xyzTrk[icoord] - xyzVtx[icoord];dst2 += di*di; }
722 t->StartTimeIntegral();
723 t->AddTimeStep(TMath::Sqrt(dst2));
725 // transfer the time integral to ESD track
726 esd->SetStatus(AliESDtrack::kTIME);
727 Double_t times[10];t->GetIntegratedTimes(times); esd->SetIntegratedTimes(times);
728 esd->SetIntegratedLength(t->GetIntegratedLength());
730 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
732 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
733 ResetTrackToFollow(*t);
735 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
736 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
737 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
741 fTrackToFollow.SetLabel(t->GetLabel());
742 //fTrackToFollow.CookdEdx();
743 CookLabel(&fTrackToFollow,0.); //For comparison only
744 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
745 //UseClusters(&fTrackToFollow);
751 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
753 fTrackingPhase="Default";
757 //------------------------------------------------------------------------
758 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
759 //--------------------------------------------------------------------
760 // This functions refits ITS tracks using the
761 // "inward propagated" TPC tracks
762 // The clusters must be loaded !
763 //--------------------------------------------------------------------
764 fTrackingPhase="RefitInward";
766 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
768 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
769 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
771 Int_t nentr=event->GetNumberOfTracks();
772 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
774 // only for PlaneEff and in case of SPD (for FO studies)
775 if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
776 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 &&
777 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) {
778 for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
782 for (Int_t i=0; i<nentr; i++) {
783 AliESDtrack *esd=event->GetTrack(i);
785 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
786 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
787 if (esd->GetStatus()&AliESDtrack::kTPCout)
788 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
790 AliITStrackMI *t = new AliITStrackMI(*esd);
792 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
793 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
798 ResetTrackToFollow(*t);
799 fTrackToFollow.ResetClusters();
801 // ITS standalone tracks
802 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
803 fTrackToFollow.ResetCovariance(10.);
804 // protection for loopers that can have parameters screwed up
805 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
806 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
813 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
814 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
816 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
817 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
818 AliDebug(2," refit OK");
819 fTrackToFollow.SetLabel(t->GetLabel());
820 // fTrackToFollow.CookdEdx();
821 CookdEdx(&fTrackToFollow);
823 CookLabel(&fTrackToFollow,0.0); //For comparison only
826 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
827 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
828 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
829 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
830 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
831 Double_t r[3]={0.,0.,0.};
833 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
840 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
842 fTrackingPhase="Default";
846 //------------------------------------------------------------------------
847 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
848 //--------------------------------------------------------------------
849 // Return pointer to a given cluster
850 //--------------------------------------------------------------------
851 Int_t l=(index & 0xf0000000) >> 28;
852 Int_t c=(index & 0x0fffffff) >> 00;
853 return fgLayers[l].GetCluster(c);
855 //------------------------------------------------------------------------
856 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
857 //--------------------------------------------------------------------
858 // Get track space point with index i
859 //--------------------------------------------------------------------
861 Int_t l=(index & 0xf0000000) >> 28;
862 Int_t c=(index & 0x0fffffff) >> 00;
863 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
864 Int_t idet = cl->GetDetectorIndex();
868 cl->GetGlobalXYZ(xyz);
869 cl->GetGlobalCov(cov);
871 p.SetCharge(cl->GetQ());
872 p.SetDriftTime(cl->GetDriftTime());
873 p.SetChargeRatio(cl->GetChargeRatio());
874 p.SetClusterType(cl->GetClusterType());
875 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
878 iLayer = AliGeomManager::kSPD1;
881 iLayer = AliGeomManager::kSPD2;
884 iLayer = AliGeomManager::kSDD1;
887 iLayer = AliGeomManager::kSDD2;
890 iLayer = AliGeomManager::kSSD1;
893 iLayer = AliGeomManager::kSSD2;
896 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
899 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
900 p.SetVolumeID((UShort_t)volid);
903 //------------------------------------------------------------------------
904 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
905 AliTrackPoint& p, const AliESDtrack *t) {
906 //--------------------------------------------------------------------
907 // Get track space point with index i
908 // (assign error estimated during the tracking)
909 //--------------------------------------------------------------------
911 Int_t l=(index & 0xf0000000) >> 28;
912 Int_t c=(index & 0x0fffffff) >> 00;
913 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
914 Int_t idet = cl->GetDetectorIndex();
916 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
918 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
920 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
921 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
922 Double_t alpha = t->GetAlpha();
923 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
924 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
925 phi += alpha-det.GetPhi();
926 Float_t tgphi = TMath::Tan(phi);
928 Float_t tgl = t->GetTgl(); // tgl about const along track
929 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
931 Float_t errtrky,errtrkz,covyz;
932 Bool_t addMisalErr=kFALSE;
933 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
937 cl->GetGlobalXYZ(xyz);
938 // cl->GetGlobalCov(cov);
939 Float_t pos[3] = {0.,0.,0.};
940 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
941 tmpcl.GetGlobalCov(cov);
944 p.SetCharge(cl->GetQ());
945 p.SetDriftTime(cl->GetDriftTime());
946 p.SetChargeRatio(cl->GetChargeRatio());
947 p.SetClusterType(cl->GetClusterType());
949 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
952 iLayer = AliGeomManager::kSPD1;
955 iLayer = AliGeomManager::kSPD2;
958 iLayer = AliGeomManager::kSDD1;
961 iLayer = AliGeomManager::kSDD2;
964 iLayer = AliGeomManager::kSSD1;
967 iLayer = AliGeomManager::kSSD2;
970 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
973 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
975 p.SetVolumeID((UShort_t)volid);
978 //------------------------------------------------------------------------
979 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
981 //--------------------------------------------------------------------
982 // Follow prolongation tree
983 //--------------------------------------------------------------------
985 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
986 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
989 AliESDtrack * esd = otrack->GetESDtrack();
990 if (esd->GetV0Index(0)>0) {
991 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
992 // mapping of ESD track is different as ITS track in Containers
993 // Need something more stable
994 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
995 for (Int_t i=0;i<3;i++){
996 Int_t index = esd->GetV0Index(i);
998 AliESDv0 * vertex = fEsd->GetV0(index);
999 if (vertex->GetStatus()<0) continue; // rejected V0
1001 if (esd->GetSign()>0) {
1002 vertex->SetIndex(0,esdindex);
1004 vertex->SetIndex(1,esdindex);
1008 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1010 bestarray = new TObjArray(5);
1011 bestarray->SetOwner();
1012 fBestHypothesys.AddAt(bestarray,esdindex);
1016 //setup tree of the prolongations
1018 const int kMaxTr = 100; //RS
1019 static AliITStrackMI tracks[7][kMaxTr];
1020 AliITStrackMI *currenttrack;
1021 static AliITStrackMI currenttrack1;
1022 static AliITStrackMI currenttrack2;
1023 static AliITStrackMI backuptrack;
1025 Int_t nindexes[7][kMaxTr];
1026 Float_t normalizedchi2[kMaxTr];
1027 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1028 otrack->SetNSkipped(0);
1029 new (&(tracks[6][0])) AliITStrackMI(*otrack);
1031 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1032 Int_t modstatus = 1; // found
1036 // follow prolongations
1037 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1038 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1041 AliITSlayer &layer=fgLayers[ilayer];
1042 Double_t r = layer.GetR();
1048 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1049 // printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1051 if (ntracks[ilayer]>=kMaxTr) break;
1052 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1053 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1054 if (ntracks[ilayer]>15+ilayer){
1055 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1056 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1059 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1061 // material between SSD and SDD, SDD and SPD
1063 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1065 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1069 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1070 Int_t idet=layer.FindDetectorIndex(phi,z);
1072 Double_t trackGlobXYZ1[3];
1073 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1075 // Get the budget to the primary vertex for the current track being prolonged
1076 Double_t budgetToPrimVertex = 0;
1077 double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1080 if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1081 else budgetToPrimVertex = GetEffectiveThickness();
1083 // check if we allow a prolongation without point
1084 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1086 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1087 // propagate to the layer radius
1088 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1089 if(!vtrack->Propagate(xToGo)) continue;
1090 // apply correction for material of the current layer
1091 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1092 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1093 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1094 vtrack->SetClIndex(ilayer,-1);
1095 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1096 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1097 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1099 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1101 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1102 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1108 // track outside layer acceptance in z
1109 if (idet<0) continue;
1111 //propagate to the intersection with the detector plane
1112 const AliITSdetector &det=layer.GetDetector(idet);
1113 new(¤ttrack2) AliITStrackMI(currenttrack1);
1114 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1115 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1116 currenttrack1.SetDetectorIndex(idet);
1117 currenttrack2.SetDetectorIndex(idet);
1118 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1121 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1123 // road in global (rphi,z) [i.e. in tracking ref. system]
1124 Double_t zmin,zmax,ymin,ymax;
1125 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1127 // select clusters in road
1128 layer.SelectClusters(zmin,zmax,ymin,ymax);
1129 //********************
1131 // Define criteria for track-cluster association
1132 Double_t msz = currenttrack1.GetSigmaZ2() +
1133 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1134 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1135 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1136 Double_t msy = currenttrack1.GetSigmaY2() +
1137 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1138 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1139 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1141 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1142 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1144 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1145 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1147 msz = 1./msz; // 1/RoadZ^2
1148 msy = 1./msy; // 1/RoadY^2
1152 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1154 const AliITSRecPoint *cl=0;
1156 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1157 Bool_t deadzoneSPD=kFALSE;
1158 currenttrack = ¤ttrack1;
1160 // check if the road contains a dead zone
1161 Bool_t noClusters = kFALSE;
1162 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1163 if (noClusters) AliDebug(2,"no clusters in road");
1164 Double_t dz=0.5*(zmax-zmin);
1165 Double_t dy=0.5*(ymax-ymin);
1166 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1167 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1168 // create a prolongation without clusters (check also if there are no clusters in the road)
1171 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1172 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1173 updatetrack->SetClIndex(ilayer,-1);
1175 modstatus = 5; // no cls in road
1176 } else if (dead==1) {
1177 modstatus = 7; // holes in z in SPD
1178 } else if (dead==2 || dead==3 || dead==4) {
1179 modstatus = 2; // dead from OCDB
1181 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1182 // apply correction for material of the current layer
1183 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1184 if (constrain) { // apply vertex constrain
1185 updatetrack->SetConstrain(constrain);
1186 Bool_t isPrim = kTRUE;
1187 if (ilayer<4) { // check that it's close to the vertex
1188 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1189 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1190 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1191 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1192 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1194 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1196 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1197 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1200 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1202 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1203 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1205 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1206 updatetrack->SetDeadZoneProbability(ilayer,1.);
1207 } else if (dead==4) { // at least a single dead channel from OCDB
1208 updatetrack->SetDeadZoneProbability(ilayer,0.);
1215 // loop over clusters in the road
1216 while ((cl=layer.GetNextCluster(clidx))!=0) {
1217 if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters
1218 Bool_t changedet =kFALSE;
1219 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1220 Int_t idetc=cl->GetDetectorIndex();
1222 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1223 // take into account misalignment (bring track to real detector plane)
1224 Double_t xTrOrig = currenttrack->GetX();
1225 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1226 // a first cut on track-cluster distance
1227 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1228 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1229 { // cluster not associated to track
1230 AliDebug(2,"not associated");
1231 // MvL: added here as well
1232 // bring track back to ideal detector plane
1233 currenttrack->Propagate(xTrOrig);
1236 // bring track back to ideal detector plane
1237 if (!currenttrack->Propagate(xTrOrig)) continue;
1238 } else { // have to move track to cluster's detector
1239 const AliITSdetector &detc=layer.GetDetector(idetc);
1240 // a first cut on track-cluster distance
1242 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1243 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1244 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1245 continue; // cluster not associated to track
1247 new (&backuptrack) AliITStrackMI(currenttrack2);
1249 currenttrack =¤ttrack2;
1250 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1251 new (currenttrack) AliITStrackMI(backuptrack);
1255 currenttrack->SetDetectorIndex(idetc);
1256 // Get again the budget to the primary vertex
1257 // for the current track being prolonged, if had to change detector
1258 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1261 // calculate track-clusters chi2
1262 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1264 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1265 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1266 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1267 if (ntracks[ilayer]>=kMaxTr) continue;
1268 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1269 updatetrack->SetClIndex(ilayer,-1);
1270 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1272 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1273 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1274 AliDebug(2,"update failed");
1277 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1278 modstatus = 1; // found
1279 } else { // virtual cluster in dead zone
1280 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1281 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1282 modstatus = 7; // holes in z in SPD
1286 Float_t xlocnewdet,zlocnewdet;
1287 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1288 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1291 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1293 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1295 // apply correction for material of the current layer
1296 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1298 if (constrain) { // apply vertex constrain
1299 updatetrack->SetConstrain(constrain);
1300 Bool_t isPrim = kTRUE;
1301 if (ilayer<4) { // check that it's close to the vertex
1302 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1303 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1304 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1305 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1306 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1308 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1310 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1311 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1313 } //apply vertex constrain
1315 } // create new hypothesis
1317 AliDebug(2,"chi2 too large");
1320 } // loop over possible prolongations
1322 // allow one prolongation without clusters
1323 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1324 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1325 // apply correction for material of the current layer
1326 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1327 vtrack->SetClIndex(ilayer,-1);
1328 modstatus = 3; // skipped
1329 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1330 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1332 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1333 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1335 vtrack->IncrementNSkipped();
1340 } // loop over tracks in layer ilayer+1
1342 //loop over track candidates for the current layer
1348 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1349 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1350 if (normalizedchi2[itrack] <
1351 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1355 if (constrain) { // constrain
1356 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1358 } else { // !constrain
1359 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1364 // sort tracks by increasing normalized chi2
1365 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1366 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1367 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1368 // if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1369 if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9);
1370 } // end loop over layers
1374 // Now select tracks to be kept
1376 Int_t max = constrain ? 20 : 5;
1378 // tracks that reach layer 0 (SPD inner)
1379 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1380 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1381 if (track.GetNumberOfClusters()<2) continue;
1382 if (!constrain && track.GetNormChi2(0) >
1383 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1386 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1389 // tracks that reach layer 1 (SPD outer)
1390 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1391 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1392 if (track.GetNumberOfClusters()<4) continue;
1393 if (!constrain && track.GetNormChi2(1) >
1394 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1395 if (constrain) track.IncrementNSkipped();
1397 track.SetD(0,track.GetD(GetX(),GetY()));
1398 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1399 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1400 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1403 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1406 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1408 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1409 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1410 if (track.GetNumberOfClusters()<3) continue;
1411 if (track.GetNormChi2(2) >
1412 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1413 track.SetD(0,track.GetD(GetX(),GetY()));
1414 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1415 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1416 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1418 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1424 // register best track of each layer - important for V0 finder
1426 for (Int_t ilayer=0;ilayer<5;ilayer++){
1427 if (ntracks[ilayer]==0) continue;
1428 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1429 if (track.GetNumberOfClusters()<1) continue;
1430 CookLabel(&track,0);
1431 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1435 // update TPC V0 information
1437 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1438 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1439 for (Int_t i=0;i<3;i++){
1440 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1441 if (index==0) break;
1442 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1443 if (vertex->GetStatus()<0) continue; // rejected V0
1445 if (otrack->GetSign()>0) {
1446 vertex->SetIndex(0,esdindex);
1449 vertex->SetIndex(1,esdindex);
1451 //find nearest layer with track info
1452 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1453 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1454 Int_t nearest = nearestold;
1455 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1456 if (ntracks[nearest]==0){
1461 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1462 if (nearestold<5&&nearest<5){
1463 Bool_t accept = track.GetNormChi2(nearest)<10;
1465 if (track.GetSign()>0) {
1466 vertex->SetParamP(track);
1467 vertex->Update(fprimvertex);
1468 //vertex->SetIndex(0,track.fESDtrack->GetID());
1469 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1471 vertex->SetParamN(track);
1472 vertex->Update(fprimvertex);
1473 //vertex->SetIndex(1,track.fESDtrack->GetID());
1474 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1476 vertex->SetStatus(vertex->GetStatus()+1);
1478 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1485 //------------------------------------------------------------------------
1486 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1488 //--------------------------------------------------------------------
1491 return fgLayers[layer];
1493 //------------------------------------------------------------------------
1494 AliITStrackerMI::AliITSlayer::AliITSlayer():
1524 //--------------------------------------------------------------------
1525 //default AliITSlayer constructor
1526 //--------------------------------------------------------------------
1527 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1528 fClusterWeight[i]=0;
1529 fClusterTracks[0][i]=-1;
1530 fClusterTracks[1][i]=-1;
1531 fClusterTracks[2][i]=-1;
1532 fClusterTracks[3][i]=-1;
1539 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1540 for (Int_t j1=0; j1<6; j1++) {
1541 fClusters5[j1][j]=0;
1542 fClusterIndex5[j1][j]=-1;
1551 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1552 for (Int_t j1=0; j1<11; j1++) {
1553 fClusters10[j1][j]=0;
1554 fClusterIndex10[j1][j]=-1;
1563 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1564 for (Int_t j1=0; j1<21; j1++) {
1565 fClusters20[j1][j]=0;
1566 fClusterIndex20[j1][j]=-1;
1574 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1579 //------------------------------------------------------------------------
1580 AliITStrackerMI::AliITSlayer::
1581 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1610 //--------------------------------------------------------------------
1611 //main AliITSlayer constructor
1612 //--------------------------------------------------------------------
1613 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1614 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1616 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1617 fClusterWeight[i]=0;
1618 fClusterTracks[0][i]=-1;
1619 fClusterTracks[1][i]=-1;
1620 fClusterTracks[2][i]=-1;
1621 fClusterTracks[3][i]=-1;
1629 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1630 for (Int_t j1=0; j1<6; j1++) {
1631 fClusters5[j1][j]=0;
1632 fClusterIndex5[j1][j]=-1;
1641 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1642 for (Int_t j1=0; j1<11; j1++) {
1643 fClusters10[j1][j]=0;
1644 fClusterIndex10[j1][j]=-1;
1653 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1654 for (Int_t j1=0; j1<21; j1++) {
1655 fClusters20[j1][j]=0;
1656 fClusterIndex20[j1][j]=-1;
1664 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1670 //------------------------------------------------------------------------
1671 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1673 fPhiOffset(layer.fPhiOffset),
1674 fNladders(layer.fNladders),
1675 fZOffset(layer.fZOffset),
1676 fNdetectors(layer.fNdetectors),
1677 fDetectors(layer.fDetectors),
1682 fClustersCs(layer.fClustersCs),
1683 fClusterIndexCs(layer.fClusterIndexCs),
1687 fCurrentSlice(layer.fCurrentSlice),
1695 fAccepted(layer.fAccepted),
1697 fMaxSigmaClY(layer.fMaxSigmaClY),
1698 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1699 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1704 //------------------------------------------------------------------------
1705 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1706 //--------------------------------------------------------------------
1707 // AliITSlayer destructor
1708 //--------------------------------------------------------------------
1709 delete [] fDetectors;
1710 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1711 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1712 fClusterWeight[i]=0;
1713 fClusterTracks[0][i]=-1;
1714 fClusterTracks[1][i]=-1;
1715 fClusterTracks[2][i]=-1;
1716 fClusterTracks[3][i]=-1;
1719 //------------------------------------------------------------------------
1720 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1721 //--------------------------------------------------------------------
1722 // This function removes loaded clusters
1723 //--------------------------------------------------------------------
1724 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1725 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1726 fClusterWeight[i]=0;
1727 fClusterTracks[0][i]=-1;
1728 fClusterTracks[1][i]=-1;
1729 fClusterTracks[2][i]=-1;
1730 fClusterTracks[3][i]=-1;
1736 //------------------------------------------------------------------------
1737 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1738 //--------------------------------------------------------------------
1739 // This function reset weights of the clusters
1740 //--------------------------------------------------------------------
1741 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1742 fClusterWeight[i]=0;
1743 fClusterTracks[0][i]=-1;
1744 fClusterTracks[1][i]=-1;
1745 fClusterTracks[2][i]=-1;
1746 fClusterTracks[3][i]=-1;
1748 for (Int_t i=0; i<fN;i++) {
1749 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1750 if (cl&&cl->IsUsed()) cl->Use();
1754 //------------------------------------------------------------------------
1755 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1756 //--------------------------------------------------------------------
1757 // This function calculates the road defined by the cluster density
1758 //--------------------------------------------------------------------
1760 for (Int_t i=0; i<fN; i++) {
1761 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1763 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1765 //------------------------------------------------------------------------
1766 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1767 //--------------------------------------------------------------------
1768 //This function adds a cluster to this layer
1769 //--------------------------------------------------------------------
1770 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1776 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1778 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1779 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1780 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1781 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1782 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1783 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1786 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1787 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1788 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1789 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1793 //------------------------------------------------------------------------
1794 void AliITStrackerMI::AliITSlayer::SortClusters()
1799 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1800 Float_t *z = new Float_t[fN];
1801 Int_t * index = new Int_t[fN];
1803 fMaxSigmaClY=0.; //AD
1804 fMaxSigmaClZ=0.; //AD
1806 for (Int_t i=0;i<fN;i++){
1807 z[i] = fClusters[i]->GetZ();
1808 // save largest errors in y and z for this layer
1809 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1810 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1812 TMath::Sort(fN,z,index,kFALSE);
1813 for (Int_t i=0;i<fN;i++){
1814 clusters[i] = fClusters[index[i]];
1817 for (Int_t i=0;i<fN;i++){
1818 fClusters[i] = clusters[i];
1819 fZ[i] = fClusters[i]->GetZ();
1820 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1821 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1822 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1832 for (Int_t i=0;i<fN;i++){
1833 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1834 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1835 fClusterIndex[i] = i;
1839 fDy5 = (fYB[1]-fYB[0])/5.;
1840 fDy10 = (fYB[1]-fYB[0])/10.;
1841 fDy20 = (fYB[1]-fYB[0])/20.;
1842 for (Int_t i=0;i<6;i++) fN5[i] =0;
1843 for (Int_t i=0;i<11;i++) fN10[i]=0;
1844 for (Int_t i=0;i<21;i++) fN20[i]=0;
1846 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;}
1847 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;}
1848 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;}
1851 for (Int_t i=0;i<fN;i++)
1852 for (Int_t irot=-1;irot<=1;irot++){
1853 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1855 for (Int_t slice=0; slice<6;slice++){
1856 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1857 fClusters5[slice][fN5[slice]] = fClusters[i];
1858 fY5[slice][fN5[slice]] = curY;
1859 fZ5[slice][fN5[slice]] = fZ[i];
1860 fClusterIndex5[slice][fN5[slice]]=i;
1865 for (Int_t slice=0; slice<11;slice++){
1866 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1867 fClusters10[slice][fN10[slice]] = fClusters[i];
1868 fY10[slice][fN10[slice]] = curY;
1869 fZ10[slice][fN10[slice]] = fZ[i];
1870 fClusterIndex10[slice][fN10[slice]]=i;
1875 for (Int_t slice=0; slice<21;slice++){
1876 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1877 fClusters20[slice][fN20[slice]] = fClusters[i];
1878 fY20[slice][fN20[slice]] = curY;
1879 fZ20[slice][fN20[slice]] = fZ[i];
1880 fClusterIndex20[slice][fN20[slice]]=i;
1887 // consistency check
1889 for (Int_t i=0;i<fN-1;i++){
1895 for (Int_t slice=0;slice<21;slice++)
1896 for (Int_t i=0;i<fN20[slice]-1;i++){
1897 if (fZ20[slice][i]>fZ20[slice][i+1]){
1904 //------------------------------------------------------------------------
1905 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1906 //--------------------------------------------------------------------
1907 // This function returns the index of the nearest cluster
1908 //--------------------------------------------------------------------
1911 if (fCurrentSlice<0) {
1920 if (ncl==0) return 0;
1921 Int_t b=0, e=ncl-1, m=(b+e)/2;
1922 for (; b<e; m=(b+e)/2) {
1923 // if (z > fClusters[m]->GetZ()) b=m+1;
1924 if (z > zcl[m]) b=m+1;
1929 //------------------------------------------------------------------------
1930 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 {
1931 //--------------------------------------------------------------------
1932 // This function computes the rectangular road for this track
1933 //--------------------------------------------------------------------
1936 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1937 // take into account the misalignment: propagate track to misaligned detector plane
1938 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1940 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1941 TMath::Sqrt(track->GetSigmaZ2() +
1942 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1943 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1944 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1945 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1946 TMath::Sqrt(track->GetSigmaY2() +
1947 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1948 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1949 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1951 // track at boundary between detectors, enlarge road
1952 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1953 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1954 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1955 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1956 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1957 Float_t tgl = TMath::Abs(track->GetTgl());
1958 if (tgl > 1.) tgl=1.;
1959 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1960 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1961 Float_t snp = TMath::Abs(track->GetSnp());
1962 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1963 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1966 // add to the road a term (up to 2-3 mm) to deal with misalignments
1967 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1968 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1970 Double_t r = fgLayers[ilayer].GetR();
1971 zmin = track->GetZ() - dz;
1972 zmax = track->GetZ() + dz;
1973 ymin = track->GetY() + r*det.GetPhi() - dy;
1974 ymax = track->GetY() + r*det.GetPhi() + dy;
1976 // bring track back to idead detector plane
1977 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1981 //------------------------------------------------------------------------
1982 void AliITStrackerMI::AliITSlayer::
1983 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1984 //--------------------------------------------------------------------
1985 // This function sets the "window"
1986 //--------------------------------------------------------------------
1988 Double_t circle=2*TMath::Pi()*fR;
1994 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1995 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1996 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1997 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1998 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
2000 Float_t ymiddle = (fYmin+fYmax)*0.5;
2001 if (ymiddle<fYB[0]) {
2002 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
2003 } else if (ymiddle>fYB[1]) {
2004 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
2010 fClustersCs = fClusters;
2011 fClusterIndexCs = fClusterIndex;
2017 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2018 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2019 if (slice<0) slice=0;
2020 if (slice>20) slice=20;
2021 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2023 fCurrentSlice=slice;
2024 fClustersCs = fClusters20[fCurrentSlice];
2025 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2026 fYcs = fY20[fCurrentSlice];
2027 fZcs = fZ20[fCurrentSlice];
2028 fNcs = fN20[fCurrentSlice];
2033 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2034 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2035 if (slice<0) slice=0;
2036 if (slice>10) slice=10;
2037 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2039 fCurrentSlice=slice;
2040 fClustersCs = fClusters10[fCurrentSlice];
2041 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2042 fYcs = fY10[fCurrentSlice];
2043 fZcs = fZ10[fCurrentSlice];
2044 fNcs = fN10[fCurrentSlice];
2049 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2050 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2051 if (slice<0) slice=0;
2052 if (slice>5) slice=5;
2053 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2055 fCurrentSlice=slice;
2056 fClustersCs = fClusters5[fCurrentSlice];
2057 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2058 fYcs = fY5[fCurrentSlice];
2059 fZcs = fZ5[fCurrentSlice];
2060 fNcs = fN5[fCurrentSlice];
2064 fI = FindClusterIndex(fZmin);
2065 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2071 //------------------------------------------------------------------------
2072 Int_t AliITStrackerMI::AliITSlayer::
2073 FindDetectorIndex(Double_t phi, Double_t z) const {
2074 //--------------------------------------------------------------------
2075 //This function finds the detector crossed by the track
2076 //--------------------------------------------------------------------
2078 if (fZOffset<0) // old geometry
2079 dphi = -(phi-fPhiOffset);
2080 else // new geometry
2081 dphi = phi-fPhiOffset;
2084 if (dphi < 0) dphi += 2*TMath::Pi();
2085 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2086 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2087 if (np>=fNladders) np-=fNladders;
2088 if (np<0) np+=fNladders;
2091 Double_t dz=fZOffset-z;
2092 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2093 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2094 if (nz>=fNdetectors || nz<0) {
2095 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2099 // ad hoc correction for 3rd ladder of SDD inner layer,
2100 // which is reversed (rotated by pi around local y)
2101 // this correction is OK only from AliITSv11Hybrid onwards
2102 if (GetR()>12. && GetR()<20.) { // SDD inner
2103 if(np==2) { // 3rd ladder
2104 Double_t posMod252[3];
2105 AliITSgeomTGeo::GetTranslation(252,posMod252);
2106 // check the Z coordinate of Mod 252: if negative
2107 // (old SDD geometry in AliITSv11Hybrid)
2108 // the swap of numeration whould be applied
2109 if(posMod252[2]<0.){
2110 nz = (fNdetectors-1) - nz;
2114 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2117 return np*fNdetectors + nz;
2119 //------------------------------------------------------------------------
2120 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2122 //--------------------------------------------------------------------
2123 // This function returns clusters within the "window"
2124 //--------------------------------------------------------------------
2126 if (fCurrentSlice<0) {
2127 Double_t rpi2 = 2.*fR*TMath::Pi();
2128 for (Int_t i=fI; i<fImax; i++) {
2131 if (fYmax<y) y -= rpi2;
2132 if (fYmin>y) y += rpi2;
2133 if (y<fYmin) continue;
2134 if (y>fYmax) continue;
2136 // skip clusters that are in "extended" road but they
2137 // 3sigma error does not touch the original road
2138 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2139 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2141 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2144 return fClusters[i];
2147 for (Int_t i=fI; i<fImax; i++) {
2148 if (fYcs[i]<fYmin) continue;
2149 if (fYcs[i]>fYmax) continue;
2150 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2151 ci=fClusterIndexCs[i];
2153 return fClustersCs[i];
2158 //------------------------------------------------------------------------
2159 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2161 //--------------------------------------------------------------------
2162 // This function returns the layer thickness at this point (units X0)
2163 //--------------------------------------------------------------------
2165 x0=AliITSRecoParam::GetX0Air();
2166 if (43<fR&&fR<45) { //SSD2
2169 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2170 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2171 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2172 for (Int_t i=0; i<12; i++) {
2173 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2174 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2178 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2179 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2183 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2184 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2187 if (37<fR&&fR<41) { //SSD1
2190 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2191 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2192 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2193 for (Int_t i=0; i<11; i++) {
2194 if (TMath::Abs(z-3.9*i)<0.15) {
2195 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2199 if (TMath::Abs(z+3.9*i)<0.15) {
2200 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2204 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2205 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2208 if (13<fR&&fR<26) { //SDD
2211 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2213 if (TMath::Abs(y-1.80)<0.55) {
2215 for (Int_t j=0; j<20; j++) {
2216 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2217 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2220 if (TMath::Abs(y+1.80)<0.55) {
2222 for (Int_t j=0; j<20; j++) {
2223 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2224 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2228 for (Int_t i=0; i<4; i++) {
2229 if (TMath::Abs(z-7.3*i)<0.60) {
2231 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2234 if (TMath::Abs(z+7.3*i)<0.60) {
2236 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2241 if (6<fR&&fR<8) { //SPD2
2242 Double_t dd=0.0063; x0=21.5;
2244 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2245 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2247 if (3<fR&&fR<5) { //SPD1
2248 Double_t dd=0.0063; x0=21.5;
2250 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2251 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2256 //------------------------------------------------------------------------
2257 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2259 fRmisal(det.fRmisal),
2261 fSinPhi(det.fSinPhi),
2262 fCosPhi(det.fCosPhi),
2268 fNChips(det.fNChips),
2269 fChipIsBad(det.fChipIsBad)
2273 //------------------------------------------------------------------------
2274 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2275 const AliITSDetTypeRec *detTypeRec)
2277 //--------------------------------------------------------------------
2278 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2279 //--------------------------------------------------------------------
2281 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2282 // while in the tracker they start from 0 for each layer
2283 for(Int_t il=0; il<ilayer; il++)
2284 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2287 if (ilayer==0 || ilayer==1) { // ---------- SPD
2289 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2291 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2294 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2298 // Get calibration from AliITSDetTypeRec
2299 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2300 calib->SetModuleIndex(idet);
2301 AliITSCalibration *calibSPDdead = 0;
2302 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2303 if (calib->IsBad() ||
2304 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2307 // printf("lay %d bad %d\n",ilayer,idet);
2310 // Get segmentation from AliITSDetTypeRec
2311 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2313 // Read info about bad chips
2314 fNChips = segm->GetMaximumChipIndex()+1;
2315 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2316 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2317 fChipIsBad = new Bool_t[fNChips];
2318 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2319 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2320 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2321 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2326 //------------------------------------------------------------------------
2327 Double_t AliITStrackerMI::GetEffectiveThickness()
2329 //--------------------------------------------------------------------
2330 // Returns the thickness between the current layer and the vertex (units X0)
2331 //--------------------------------------------------------------------
2334 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2335 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2336 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2340 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2341 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2345 Double_t xn=fgLayers[fI].GetR();
2346 for (Int_t i=0; i<fI; i++) {
2347 Double_t xi=fgLayers[i].GetR();
2348 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2354 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2355 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2358 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2359 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2364 //------------------------------------------------------------------------
2365 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2367 //--------------------------------------------------------------------
2368 // Returns the array of layers between the current layer and the vertex
2369 //--------------------------------------------------------------------
2372 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2373 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2374 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2379 for (int il=fI;il--;) {
2382 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2383 xMS[nl++] = AliITSRecoParam::GetrInsideShield(1);
2386 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2387 xMS[nl++] = AliITSRecoParam::GetrInsideShield(0);
2390 x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2391 xMS[nl++] = fgLayers[il].GetR();
2396 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2397 xMS[nl++] = AliITSRecoParam::GetrPipe();
2403 //------------------------------------------------------------------------
2404 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2405 //-------------------------------------------------------------------
2406 // This function returns number of clusters within the "window"
2407 //--------------------------------------------------------------------
2409 for (Int_t i=fI; i<fN; i++) {
2410 const AliITSRecPoint *c=fClusters[i];
2411 if (c->GetZ() > fZmax) break;
2412 if (c->IsUsed()) continue;
2413 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2414 Double_t y=fR*det.GetPhi() + c->GetY();
2416 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2417 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2419 if (y<fYmin) continue;
2420 if (y>fYmax) continue;
2425 //------------------------------------------------------------------------
2426 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2427 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2429 //--------------------------------------------------------------------
2430 // This function refits the track "track" at the position "x" using
2431 // the clusters from "clusters"
2432 // If "extra"==kTRUE,
2433 // the clusters from overlapped modules get attached to "track"
2434 // If "planeff"==kTRUE,
2435 // special approach for plane efficiency evaluation is applyed
2436 //--------------------------------------------------------------------
2438 Int_t index[AliITSgeomTGeo::kNLayers];
2440 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2441 Int_t nc=clusters->GetNumberOfClusters();
2442 for (k=0; k<nc; k++) {
2443 Int_t idx=clusters->GetClusterIndex(k);
2444 Int_t ilayer=(idx&0xf0000000)>>28;
2448 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2450 //------------------------------------------------------------------------
2451 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2452 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2454 //--------------------------------------------------------------------
2455 // This function refits the track "track" at the position "x" using
2456 // the clusters from array
2457 // If "extra"==kTRUE,
2458 // the clusters from overlapped modules get attached to "track"
2459 // If "planeff"==kTRUE,
2460 // special approach for plane efficiency evaluation is applyed
2461 //--------------------------------------------------------------------
2462 Int_t index[AliITSgeomTGeo::kNLayers];
2464 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2466 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2467 index[k]=clusters[k];
2470 // special for cosmics and TPC prolonged tracks:
2471 // propagate to the innermost of:
2472 // - innermost layer crossed by the track
2473 // - innermost layer where a cluster was associated to the track
2474 static AliITSRecoParam *repa = NULL;
2476 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2478 repa = AliITSRecoParam::GetHighFluxParam();
2479 AliWarning("Using default AliITSRecoParam class");
2482 Int_t evsp=repa->GetEventSpecie();
2484 if(track->GetESDtrack()) trStatus=track->GetStatus();
2485 Int_t innermostlayer=0;
2486 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2488 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2489 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2490 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2491 index[innermostlayer] >= 0 ) break;
2494 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2497 Int_t modstatus=1; // found
2499 Int_t from, to, step;
2500 if (xx > track->GetX()) {
2501 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2504 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2507 TString dir = (step>0 ? "outward" : "inward");
2509 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2510 AliITSlayer &layer=fgLayers[ilayer];
2511 Double_t r=layer.GetR();
2513 if (step<0 && xx>r) break;
2515 // material between SSD and SDD, SDD and SPD
2516 Double_t hI=ilayer-0.5*step;
2517 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2518 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2519 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2520 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2523 Double_t oldGlobXYZ[3];
2524 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2526 // continue if we are already beyond this layer
2527 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2528 if(step>0 && oldGlobR > r) continue; // going outward
2529 if(step<0 && oldGlobR < r) continue; // going inward
2532 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2534 Int_t idet=layer.FindDetectorIndex(phi,z);
2536 // check if we allow a prolongation without point for large-eta tracks
2537 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2539 modstatus = 4; // out in z
2540 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2541 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2544 // apply correction for material of the current layer
2545 // add time if going outward
2546 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2550 if (idet<0) return kFALSE;
2553 const AliITSdetector &det=layer.GetDetector(idet);
2554 // only for ITS-SA tracks refit
2555 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2557 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2559 track->SetDetectorIndex(idet);
2560 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2562 Double_t dz,zmin,zmax,dy,ymin,ymax;
2564 const AliITSRecPoint *clAcc=0;
2565 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2567 Int_t idx=index[ilayer];
2568 if (idx>=0) { // cluster in this layer
2569 modstatus = 6; // no refit
2570 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2572 if (idet != cl->GetDetectorIndex()) {
2573 idet=cl->GetDetectorIndex();
2574 const AliITSdetector &detc=layer.GetDetector(idet);
2575 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2576 track->SetDetectorIndex(idet);
2577 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2579 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2580 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2584 modstatus = 1; // found
2589 } else { // no cluster in this layer
2591 modstatus = 3; // skipped
2592 // Plane Eff determination:
2593 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2594 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2595 UseTrackForPlaneEff(track,ilayer);
2598 modstatus = 5; // no cls in road
2600 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2601 dz = 0.5*(zmax-zmin);
2602 dy = 0.5*(ymax-ymin);
2603 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2604 if (dead==1) modstatus = 7; // holes in z in SPD
2605 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2610 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2611 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2613 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2616 if (extra && clAcc) { // search for extra clusters in overlapped modules
2617 AliITStrackV2 tmp(*track);
2618 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2619 layer.SelectClusters(zmin,zmax,ymin,ymax);
2621 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2623 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2624 Double_t tolerance=0.1;
2625 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2626 // only clusters in another module! (overlaps)
2627 idetExtra = clExtra->GetDetectorIndex();
2628 if (idet == idetExtra) continue;
2630 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2632 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2633 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2634 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2635 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2637 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2638 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2641 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2642 track->SetExtraModule(ilayer,idetExtra);
2644 } // end search for extra clusters in overlapped modules
2646 // Correct for material of the current layer
2648 // add time if going outward
2649 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2650 track->SetCheckInvariant(kTRUE);
2651 } // end loop on layers
2653 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2657 //------------------------------------------------------------------------
2658 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2661 // calculate normalized chi2
2662 // return NormalizedChi2(track,0);
2665 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2666 // track->fdEdxMismatch=0;
2667 Float_t dedxmismatch =0;
2668 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2670 for (Int_t i = 0;i<6;i++){
2671 if (track->GetClIndex(i)>=0){
2672 Float_t cerry, cerrz;
2673 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2675 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2678 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2679 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2680 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2682 cchi2+=(0.5-ratio)*10.;
2683 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2684 dedxmismatch+=(0.5-ratio)*10.;
2688 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2689 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2690 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2691 if (i<2) chi2+=2*cl->GetDeltaProbability();
2697 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2698 track->SetdEdxMismatch(dedxmismatch);
2702 for (Int_t i = 0;i<4;i++){
2703 if (track->GetClIndex(i)>=0){
2704 Float_t cerry, cerrz;
2705 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2706 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2709 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2710 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2714 for (Int_t i = 4;i<6;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 Float_t cerryb, cerrzb;
2722 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2723 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2726 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2727 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2732 if (track->GetESDtrack()->GetTPCsignal()>85){
2733 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2735 chi2+=(0.5-ratio)*5.;
2738 chi2+=(ratio-2.0)*3;
2742 Double_t match = TMath::Sqrt(track->GetChi22());
2743 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2744 if (!track->GetConstrain()) {
2745 if (track->GetNumberOfClusters()>2) {
2746 match/=track->GetNumberOfClusters()-2.;
2751 if (match<0) match=0;
2753 // penalty factor for missing points (NDeadZone>0), but no penalty
2754 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2755 // or there is a dead from OCDB)
2756 Float_t deadzonefactor = 0.;
2757 if (track->GetNDeadZone()>0.) {
2758 Int_t sumDeadZoneProbability=0;
2759 for(Int_t ilay=0;ilay<6;ilay++) {
2760 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2762 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2763 if(nDeadZoneWithProbNot1>0) {
2764 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2765 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2766 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2768 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2769 deadzonefactor = 3.*(1.1-deadZoneProbability);
2773 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2774 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2775 1./(1.+track->GetNSkipped()));
2776 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2777 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2780 //------------------------------------------------------------------------
2781 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2784 // return matching chi2 between two tracks
2785 Double_t largeChi2=1000.;
2787 AliITStrackMI track3(*track2);
2788 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2790 vec(0,0)=track1->GetY() - track3.GetY();
2791 vec(1,0)=track1->GetZ() - track3.GetZ();
2792 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2793 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2794 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2797 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2798 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2799 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2800 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2801 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2803 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2804 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2805 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2806 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2808 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2809 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2810 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2812 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2813 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2815 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2818 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2819 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2822 //------------------------------------------------------------------------
2823 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2826 // return probability that given point (characterized by z position and error)
2827 // is in SPD dead zone
2828 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2830 Double_t probability = 0.;
2831 Double_t nearestz = 0.,distz=0.;
2832 Int_t nearestzone = -1;
2833 Double_t mindistz = 1000.;
2835 // find closest dead zone
2836 for (Int_t i=0; i<3; i++) {
2837 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2838 if (distz<mindistz) {
2840 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2845 // too far from dead zone
2846 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2849 Double_t zmin, zmax;
2850 if (nearestzone==0) { // dead zone at z = -7
2851 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2852 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2853 } else if (nearestzone==1) { // dead zone at z = 0
2854 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2855 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2856 } else if (nearestzone==2) { // dead zone at z = +7
2857 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2858 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2863 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2865 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2866 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2867 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2870 //------------------------------------------------------------------------
2871 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2874 // calculate normalized chi2
2876 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2878 for (Int_t i = 0;i<6;i++){
2879 if (TMath::Abs(track->GetDy(i))>0){
2880 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2881 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2884 else{chi2[i]=10000;}
2887 TMath::Sort(6,chi2,index,kFALSE);
2888 Float_t max = float(ncl)*fac-1.;
2889 Float_t sumchi=0, sumweight=0;
2890 for (Int_t i=0;i<max+1;i++){
2891 Float_t weight = (i<max)?1.:(max+1.-i);
2892 sumchi+=weight*chi2[index[i]];
2895 Double_t normchi2 = sumchi/sumweight;
2898 //------------------------------------------------------------------------
2899 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2902 // calculate normalized chi2
2903 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2906 for (Int_t i=0;i<6;i++){
2907 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2908 Double_t sy1 = forwardtrack->GetSigmaY(i);
2909 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2910 Double_t sy2 = backtrack->GetSigmaY(i);
2911 Double_t sz2 = backtrack->GetSigmaZ(i);
2912 if (i<2){ sy2=1000.;sz2=1000;}
2914 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2915 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2917 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2918 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2920 res+= nz0*nz0+ny0*ny0;
2923 if (npoints>1) return
2924 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2925 //2*forwardtrack->fNUsed+
2926 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2927 1./(1.+forwardtrack->GetNSkipped()));
2930 //------------------------------------------------------------------------
2931 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2932 //--------------------------------------------------------------------
2933 // Return pointer to a given cluster
2934 //--------------------------------------------------------------------
2935 Int_t l=(index & 0xf0000000) >> 28;
2936 Int_t c=(index & 0x0fffffff) >> 00;
2937 return fgLayers[l].GetWeight(c);
2939 //------------------------------------------------------------------------
2940 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2942 //---------------------------------------------
2943 // register track to the list
2945 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2948 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2949 Int_t index = track->GetClusterIndex(icluster);
2950 Int_t l=(index & 0xf0000000) >> 28;
2951 Int_t c=(index & 0x0fffffff) >> 00;
2952 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2953 for (Int_t itrack=0;itrack<4;itrack++){
2954 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2955 fgLayers[l].SetClusterTracks(itrack,c,id);
2961 //------------------------------------------------------------------------
2962 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2964 //---------------------------------------------
2965 // unregister track from the list
2966 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2967 Int_t index = track->GetClusterIndex(icluster);
2968 Int_t l=(index & 0xf0000000) >> 28;
2969 Int_t c=(index & 0x0fffffff) >> 00;
2970 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2971 for (Int_t itrack=0;itrack<4;itrack++){
2972 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2973 fgLayers[l].SetClusterTracks(itrack,c,-1);
2978 //------------------------------------------------------------------------
2979 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2981 //-------------------------------------------------------------
2982 //get number of shared clusters
2983 //-------------------------------------------------------------
2985 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2986 // mean number of clusters
2987 Float_t *ny = GetNy(id), *nz = GetNz(id);
2990 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2991 Int_t index = track->GetClusterIndex(icluster);
2992 Int_t l=(index & 0xf0000000) >> 28;
2993 Int_t c=(index & 0x0fffffff) >> 00;
2994 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2995 // if (ny[l]<1.e-13){
2996 // printf("problem\n");
2998 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3002 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3003 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3004 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3006 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3009 deltan = (cl->GetNz()-nz[l]);
3011 if (deltan>2.0) continue; // extended - highly probable shared cluster
3012 weight = 2./TMath::Max(3.+deltan,2.);
3014 for (Int_t itrack=0;itrack<4;itrack++){
3015 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3017 clist[l] = (AliITSRecPoint*)GetCluster(index);
3018 track->SetSharedWeight(l,weight);
3024 track->SetNUsed(shared);
3027 //------------------------------------------------------------------------
3028 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3031 // find first shared track
3033 // mean number of clusters
3034 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
3036 for (Int_t i=0;i<6;i++) overlist[i]=-1;
3037 Int_t sharedtrack=100000;
3038 Int_t tracks[24],trackindex=0;
3039 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3041 for (Int_t icluster=0;icluster<6;icluster++){
3042 if (clusterlist[icluster]<0) continue;
3043 Int_t index = clusterlist[icluster];
3044 Int_t l=(index & 0xf0000000) >> 28;
3045 Int_t c=(index & 0x0fffffff) >> 00;
3046 // if (ny[l]<1.e-13){
3047 // printf("problem\n");
3049 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3050 //if (l>3) continue;
3051 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3054 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3055 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3056 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3058 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3061 deltan = (cl->GetNz()-nz[l]);
3063 if (deltan>2.0) continue; // extended - highly probable shared cluster
3065 for (Int_t itrack=3;itrack>=0;itrack--){
3066 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3067 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3068 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
3073 if (trackindex==0) return -1;
3075 sharedtrack = tracks[0];
3077 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3080 Int_t tracks2[24], cluster[24];
3081 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
3084 for (Int_t i=0;i<trackindex;i++){
3085 if (tracks[i]<0) continue;
3086 tracks2[index] = tracks[i];
3088 for (Int_t j=i+1;j<trackindex;j++){
3089 if (tracks[j]<0) continue;
3090 if (tracks[j]==tracks[i]){
3098 for (Int_t i=0;i<index;i++){
3099 if (cluster[index]>max) {
3100 sharedtrack=tracks2[index];
3107 if (sharedtrack>=100000) return -1;
3109 // make list of overlaps
3111 for (Int_t icluster=0;icluster<6;icluster++){
3112 if (clusterlist[icluster]<0) continue;
3113 Int_t index = clusterlist[icluster];
3114 Int_t l=(index & 0xf0000000) >> 28;
3115 Int_t c=(index & 0x0fffffff) >> 00;
3116 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3117 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3119 if (cl->GetNy()>2) continue;
3120 if (cl->GetNz()>2) continue;
3123 if (cl->GetNy()>3) continue;
3124 if (cl->GetNz()>3) continue;
3127 for (Int_t itrack=3;itrack>=0;itrack--){
3128 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3129 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3137 //------------------------------------------------------------------------
3138 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3140 // try to find track hypothesys without conflicts
3141 // with minimal chi2;
3142 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3143 Int_t entries1 = arr1->GetEntriesFast();
3144 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3145 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3146 Int_t entries2 = arr2->GetEntriesFast();
3147 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3149 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3150 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3151 if (track10->Pt()>0.5+track20->Pt()) return track10;
3153 for (Int_t itrack=0;itrack<entries1;itrack++){
3154 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3155 UnRegisterClusterTracks(track,trackID1);
3158 for (Int_t itrack=0;itrack<entries2;itrack++){
3159 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3160 UnRegisterClusterTracks(track,trackID2);
3164 Float_t maxconflicts=6;
3165 Double_t maxchi2 =1000.;
3167 // get weight of hypothesys - using best hypothesy
3170 Int_t list1[6],list2[6];
3171 AliITSRecPoint *clist1[6], *clist2[6] ;
3172 RegisterClusterTracks(track10,trackID1);
3173 RegisterClusterTracks(track20,trackID2);
3174 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3175 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3176 UnRegisterClusterTracks(track10,trackID1);
3177 UnRegisterClusterTracks(track20,trackID2);
3180 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3181 Float_t nerry[6],nerrz[6];
3182 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3183 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3184 for (Int_t i=0;i<6;i++){
3185 if ( (erry1[i]>0) && (erry2[i]>0)) {
3186 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3187 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3189 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3190 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3192 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3193 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3194 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3197 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3198 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3199 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3207 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3208 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3209 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3210 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3212 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3213 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3214 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3216 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3217 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3218 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3221 Double_t sumw = w1+w2;
3225 w1 = (d2+0.5)/(d1+d2+1.);
3226 w2 = (d1+0.5)/(d1+d2+1.);
3228 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3229 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3231 // get pair of "best" hypothesys
3233 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3234 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3236 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3237 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3238 //if (track1->fFakeRatio>0) continue;
3239 RegisterClusterTracks(track1,trackID1);
3240 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3241 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3243 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3244 //if (track2->fFakeRatio>0) continue;
3246 RegisterClusterTracks(track2,trackID2);
3247 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3248 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3249 UnRegisterClusterTracks(track2,trackID2);
3251 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3252 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3253 if (nskipped>0.5) continue;
3255 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3256 if (conflict1+1<cconflict1) continue;
3257 if (conflict2+1<cconflict2) continue;
3261 for (Int_t i=0;i<6;i++){
3263 Float_t c1 =0.; // conflict coeficients
3265 if (clist1[i]&&clist2[i]){
3268 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3271 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3273 c1 = 2./TMath::Max(3.+deltan,2.);
3274 c2 = 2./TMath::Max(3.+deltan,2.);
3280 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3283 deltan = (clist1[i]->GetNz()-nz1[i]);
3285 c1 = 2./TMath::Max(3.+deltan,2.);
3292 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3295 deltan = (clist2[i]->GetNz()-nz2[i]);
3297 c2 = 2./TMath::Max(3.+deltan,2.);
3303 if (TMath::Abs(track1->GetDy(i))>0.) {
3304 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3305 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3306 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3307 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3309 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3312 if (TMath::Abs(track2->GetDy(i))>0.) {
3313 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3314 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3315 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3316 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3319 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3321 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3322 if (chi21>0) sum+=w1;
3323 if (chi22>0) sum+=w2;
3326 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3327 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3328 Double_t normchi2 = 2*conflict+sumchi2/sum;
3329 if ( normchi2 <maxchi2 ){
3332 maxconflicts = conflict;
3336 UnRegisterClusterTracks(track1,trackID1);
3339 // if (maxconflicts<4 && maxchi2<th0){
3340 if (maxchi2<th0*2.){
3341 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3342 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3343 track1->SetChi2MIP(5,maxconflicts);
3344 track1->SetChi2MIP(6,maxchi2);
3345 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3346 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3347 track1->SetChi2MIP(8,index1);
3348 fBestTrackIndex[trackID1] =index1;
3349 UpdateESDtrack(track1, AliESDtrack::kITSin);
3350 original->SetWinner(track1);
3352 else if (track10->GetChi2MIP(0)<th1){
3353 track10->SetChi2MIP(5,maxconflicts);
3354 track10->SetChi2MIP(6,maxchi2);
3355 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3356 UpdateESDtrack(track10,AliESDtrack::kITSin);
3357 original->SetWinner(track10);
3360 for (Int_t itrack=0;itrack<entries1;itrack++){
3361 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3362 UnRegisterClusterTracks(track,trackID1);
3365 for (Int_t itrack=0;itrack<entries2;itrack++){
3366 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3367 UnRegisterClusterTracks(track,trackID2);
3370 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3371 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3372 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3373 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3374 RegisterClusterTracks(track10,trackID1);
3376 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3377 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3378 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3379 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3380 RegisterClusterTracks(track20,trackID2);
3385 //------------------------------------------------------------------------
3386 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3387 //--------------------------------------------------------------------
3388 // This function marks clusters assigned to the track
3389 //--------------------------------------------------------------------
3390 AliTracker::UseClusters(t,from);
3392 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3393 //if (c->GetQ()>2) c->Use();
3394 if (c->GetSigmaZ2()>0.1) c->Use();
3395 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3396 //if (c->GetQ()>2) c->Use();
3397 if (c->GetSigmaZ2()>0.1) c->Use();
3400 //------------------------------------------------------------------------
3401 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3403 //------------------------------------------------------------------
3404 // add track to the list of hypothesys
3405 //------------------------------------------------------------------
3407 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3408 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3410 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3412 array = new TObjArray(10);
3413 fTrackHypothesys.AddAt(array,esdindex);
3415 array->AddLast(track);
3417 //------------------------------------------------------------------------
3418 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3420 //-------------------------------------------------------------------
3421 // compress array of track hypothesys
3422 // keep only maxsize best hypothesys
3423 //-------------------------------------------------------------------
3424 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3425 if (! (fTrackHypothesys.At(esdindex)) ) return;
3426 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3427 Int_t entries = array->GetEntriesFast();
3429 //- find preliminary besttrack as a reference
3430 Float_t minchi2=10000;
3432 AliITStrackMI * besttrack=0;
3434 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3435 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3436 if (!track) continue;
3437 Float_t chi2 = NormalizedChi2(track,0);
3439 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3440 track->SetLabel(tpcLabel);
3442 track->SetFakeRatio(1.);
3443 CookLabel(track,0.); //For comparison only
3445 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3446 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3447 if (track->GetNumberOfClusters()<maxn) continue;
3448 maxn = track->GetNumberOfClusters();
3449 // if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3456 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3457 delete array->RemoveAt(itrack);
3461 if (!besttrack) return;
3464 //take errors of best track as a reference
3465 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3466 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3467 for (Int_t j=0;j<6;j++) {
3468 if (besttrack->GetClIndex(j)>=0){
3469 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3470 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3471 ny[j] = besttrack->GetNy(j);
3472 nz[j] = besttrack->GetNz(j);
3476 // calculate normalized chi2
3478 Float_t * chi2 = new Float_t[entries];
3479 Int_t * index = new Int_t[entries];
3480 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3481 for (Int_t itrack=0;itrack<entries;itrack++){
3482 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3484 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3485 double chi2t = GetNormalizedChi2(track, mode);
3486 track->SetChi2MIP(0,chi2t);
3487 if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3488 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3489 chi2[itrack] = chi2t;
3492 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3493 delete array->RemoveAt(itrack);
3499 TMath::Sort(entries,chi2,index,kFALSE);
3500 besttrack = (AliITStrackMI*)array->At(index[0]);
3501 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3502 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3503 for (Int_t j=0;j<6;j++){
3504 if (besttrack->GetClIndex(j)>=0){
3505 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3506 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3507 ny[j] = besttrack->GetNy(j);
3508 nz[j] = besttrack->GetNz(j);
3513 // calculate one more time with updated normalized errors
3514 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3515 for (Int_t itrack=0;itrack<entries;itrack++){
3516 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3518 double chi2t = GetNormalizedChi2(track, mode);
3519 track->SetChi2MIP(0,chi2t);
3520 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3521 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3522 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3523 chi2[itrack] = chi2t; //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3526 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3527 delete array->RemoveAt(itrack);
3532 entries = array->GetEntriesFast();
3536 TObjArray * newarray = new TObjArray();
3537 TMath::Sort(entries,chi2,index,kFALSE);
3538 besttrack = (AliITStrackMI*)array->At(index[0]);
3540 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3542 for (Int_t j=0;j<6;j++){
3543 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3544 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3545 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3546 ny[j] = besttrack->GetNy(j);
3547 nz[j] = besttrack->GetNz(j);
3550 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3551 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3552 Float_t minn = besttrack->GetNumberOfClusters()-3;
3554 for (Int_t i=0;i<entries;i++){
3555 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3556 if (!track) continue;
3557 if (accepted>maxcut) break;
3558 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3559 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3560 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3561 delete array->RemoveAt(index[i]);
3565 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3566 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3567 if (!shortbest) accepted++;
3569 newarray->AddLast(array->RemoveAt(index[i]));
3570 for (Int_t j=0;j<6;j++){
3572 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3573 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3574 ny[j] = track->GetNy(j);
3575 nz[j] = track->GetNz(j);
3580 delete array->RemoveAt(index[i]);
3584 delete fTrackHypothesys.RemoveAt(esdindex);
3585 fTrackHypothesys.AddAt(newarray,esdindex);
3589 delete fTrackHypothesys.RemoveAt(esdindex);
3595 //------------------------------------------------------------------------
3596 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3598 //-------------------------------------------------------------
3599 // try to find best hypothesy
3600 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3601 // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3602 //-------------------------------------------------------------
3603 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3604 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3605 if (!array) return 0;
3606 Int_t entries = array->GetEntriesFast();
3607 if (!entries) return 0;
3608 Float_t minchi2 = 100000;
3609 AliITStrackMI * besttrack=0;
3611 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3612 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3613 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3614 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3616 for (Int_t i=0;i<entries;i++){
3617 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3618 if (!track) continue;
3619 Float_t sigmarfi,sigmaz;
3620 GetDCASigma(track,sigmarfi,sigmaz);
3621 track->SetDnorm(0,sigmarfi);
3622 track->SetDnorm(1,sigmaz);
3624 track->SetChi2MIP(1,1000000);
3625 track->SetChi2MIP(2,1000000);
3626 track->SetChi2MIP(3,1000000);
3629 backtrack = new(backtrack) AliITStrackMI(*track);
3630 if (track->GetConstrain()) {
3631 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3632 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3633 if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3634 else {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3636 backtrack->ResetCovariance(10.);
3638 backtrack->ResetCovariance(10.);
3640 backtrack->ResetClusters();
3642 Double_t x = original->GetX();
3643 if (!RefitAt(x,backtrack,track)) continue;
3645 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3646 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3647 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3648 track->SetChi22(GetMatchingChi2(backtrack,original));
3650 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3651 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3652 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3655 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3657 //forward track - without constraint
3658 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3659 forwardtrack->ResetClusters();
3661 if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue; // w/o fwd track MIP03 is meaningless
3662 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3663 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3664 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3666 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3667 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3668 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3669 forwardtrack->SetD(0,track->GetD(0));
3670 forwardtrack->SetD(1,track->GetD(1));
3673 AliITSRecPoint* clist[6];
3674 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3675 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3678 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3679 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3680 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3681 track->SetChi2MIP(3,1000);
3684 Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed(); //RS
3685 if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3686 else chi2 += track->GetNUsed();
3688 for (Int_t ichi=0;ichi<5;ichi++){
3689 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3691 if (chi2 < minchi2){
3692 //besttrack = new AliITStrackMI(*forwardtrack);
3694 besttrack->SetLabel(track->GetLabel());
3695 besttrack->SetFakeRatio(track->GetFakeRatio());
3697 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3698 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3699 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3703 delete forwardtrack;
3705 if (!besttrack) return 0;
3708 for (Int_t i=0;i<entries;i++){
3709 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3711 if (!track) continue;
3713 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3714 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3715 // RS: don't apply this cut when fSelectBestMIP03 is on
3716 || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3718 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3719 delete array->RemoveAt(i);
3729 SortTrackHypothesys(esdindex,checkmax,1);
3731 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3732 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3733 besttrack = (AliITStrackMI*)array->At(0);
3734 if (!besttrack) return 0;
3735 besttrack->SetChi2MIP(8,0);
3736 fBestTrackIndex[esdindex]=0;
3737 entries = array->GetEntriesFast();
3738 AliITStrackMI *longtrack =0;
3740 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3741 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3742 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3743 if (!track->GetConstrain()) continue;
3744 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3745 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3746 if (track->GetChi2MIP(0)>4.) continue;
3747 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3750 //if (longtrack) besttrack=longtrack;
3752 // RS do shared cluster analysis here only if the new sharing analysis is not requested
3753 //RRR if (fFlagFakes) return besttrack;
3756 AliITSRecPoint * clist[6];
3757 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3758 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3759 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3760 RegisterClusterTracks(besttrack,esdindex);
3767 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3768 if (sharedtrack>=0){
3770 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3772 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3778 if (shared>2.5) return 0;
3779 if (shared>1.0) return besttrack;
3781 // Don't sign clusters if not gold track
3783 if (!besttrack->IsGoldPrimary()) return besttrack;
3784 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3786 if (fConstraint[fPass]){
3790 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3791 for (Int_t i=0;i<6;i++){
3792 Int_t index = besttrack->GetClIndex(i);
3793 if (index<0) continue;
3794 Int_t ilayer = (index & 0xf0000000) >> 28;
3795 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3796 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3798 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3799 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3800 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3801 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3802 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3803 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3805 Bool_t cansign = kTRUE;
3806 for (Int_t itrack=0;itrack<entries; itrack++){
3807 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3808 if (!track) continue;
3809 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3810 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3816 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3817 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3818 if (!c->IsUsed()) c->Use();
3824 //------------------------------------------------------------------------
3825 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3828 // get "best" hypothesys
3831 Int_t nentries = itsTracks.GetEntriesFast();
3832 for (Int_t i=0;i<nentries;i++){
3833 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3834 if (!track) continue;
3835 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3836 if (!array) continue;
3837 if (array->GetEntriesFast()<=0) continue;
3839 AliITStrackMI* longtrack=0;
3841 Float_t maxchi2=1000;
3842 for (Int_t j=0;j<array->GetEntriesFast();j++){
3843 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3844 if (!trackHyp) continue;
3845 if (trackHyp->GetGoldV0()) {
3846 longtrack = trackHyp; //gold V0 track taken
3849 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3850 Float_t chi2 = trackHyp->GetChi2MIP(0);
3851 if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3852 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3854 if (fAfterV0){ // ??? RS
3855 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3857 if (chi2 > maxchi2) continue;
3858 minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3859 if (fSelectBestMIP03) minn++; // allow next to longest to win
3866 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3867 if (!longtrack) {longtrack = besttrack;}
3868 else besttrack= longtrack;
3872 AliITSRecPoint * clist[6];
3873 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3875 track->SetNUsed(shared);
3876 track->SetNSkipped(besttrack->GetNSkipped());
3877 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3879 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3883 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3884 //if (sharedtrack==-1) sharedtrack=0;
3885 if (sharedtrack>=0) {
3886 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3889 if (besttrack&&fAfterV0) {
3890 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3891 track->SetWinner(besttrack);
3894 if (fConstraint[fPass]) {
3895 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3896 track->SetWinner(besttrack);
3898 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3899 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3900 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3907 //------------------------------------------------------------------------
3908 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3911 // RS: flag those tracks which are suxpected to have fake clusters
3913 const double kThreshPt = 0.5;
3914 AliRefArray *refArr[6];
3916 for (int i=0;i<6;i++) {
3917 int ncl = fgLayers[i].GetNumberOfClusters();
3918 refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3920 Int_t nentries = itsTracks.GetEntriesFast();
3922 // fill cluster->track associations
3923 for (Int_t itr=0;itr<nentries;itr++){
3924 AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3925 if (!track) continue;
3926 AliITStrackMI* trackITS = track->GetWinner();
3927 if (!trackITS) continue;
3928 for (int il=trackITS->GetNumberOfClusters();il--;) {
3929 int idx = trackITS->GetClusterIndex(il);
3930 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3931 // if (c>fgLayers[l].GetNumberOfClusters()) continue;
3932 refArr[l]->AddReference(c, itr);
3936 const UInt_t kMaxRef = 100;
3937 UInt_t crefs[kMaxRef];
3939 // process tracks with shared clusters
3940 for (int itr=0;itr<nentries;itr++){
3941 AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3942 AliITStrackMI* trackH0 = track0->GetWinner();
3943 if (!trackH0) continue;
3944 AliESDtrack* esd0 = track0->GetESDtrack();
3946 for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3947 int idx = trackH0->GetClusterIndex(il);
3948 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3949 ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3950 if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3951 esd0->SetITSSharedFlag(l);
3952 for (int ir=ncrefs;ir--;) {
3953 if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3954 AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3955 AliITStrackMI* trackH1 = track1->GetWinner();
3956 AliESDtrack* esd1 = track1->GetESDtrack();
3957 esd1->SetITSSharedFlag(l);
3959 double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3960 if (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3961 else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3963 // select the one with smallest chi2's product
3964 res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3965 res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3967 if (res<0) esd1->SetITSFakeFlag(); // esd0 is winner
3968 else esd0->SetITSFakeFlag(); // esd1 is winner
3975 for (int i=6;i--;) delete refArr[i];
3978 //------------------------------------------------------------------------
3979 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3980 //--------------------------------------------------------------------
3981 //This function "cooks" a track label. If label<0, this track is fake.
3982 //--------------------------------------------------------------------
3985 if (track->GetESDtrack()){
3986 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3987 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3988 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3990 track->SetChi2MIP(9,0);
3992 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3993 Int_t cindex = track->GetClusterIndex(i);
3994 Int_t l=(cindex & 0xf0000000) >> 28;
3995 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3997 for (Int_t ind=0;ind<3;ind++){
3998 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3999 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
4001 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4004 Int_t nclusters = track->GetNumberOfClusters();
4005 if (nclusters > 0) //PH Some tracks don't have any cluster
4006 track->SetFakeRatio(double(nwrong)/double(nclusters));
4007 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4008 track->SetLabel(-tpcLabel);
4010 track->SetLabel(tpcLabel);
4012 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4015 //------------------------------------------------------------------------
4016 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4018 // Fill the dE/dx in this track
4020 track->SetChi2MIP(9,0);
4021 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4022 Int_t cindex = track->GetClusterIndex(i);
4023 Int_t l=(cindex & 0xf0000000) >> 28;
4024 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4025 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4027 for (Int_t ind=0;ind<3;ind++){
4028 if (cl->GetLabel(ind)==lab) isWrong=0;
4030 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4034 track->CookdEdx(low,up);
4036 //------------------------------------------------------------------------
4037 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4039 // Create some arrays
4041 if (fCoefficients) delete []fCoefficients;
4042 fCoefficients = new Float_t[ntracks*48];
4043 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4045 //------------------------------------------------------------------------
4046 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4049 // Compute predicted chi2
4051 // Take into account the mis-alignment (bring track to cluster plane)
4052 Double_t xTrOrig=track->GetX();
4053 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4054 Float_t erry,errz,covyz;
4055 Float_t theta = track->GetTgl();
4056 Float_t phi = track->GetSnp();
4057 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4058 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4059 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()));
4060 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()));
4061 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4062 // Bring the track back to detector plane in ideal geometry
4063 // [mis-alignment will be accounted for in UpdateMI()]
4064 if (!track->Propagate(xTrOrig)) return 1000.;
4066 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4067 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4069 chi2+=0.5*TMath::Min(delta/2,2.);
4070 chi2+=2.*cluster->GetDeltaProbability();
4073 track->SetNy(layer,ny);
4074 track->SetNz(layer,nz);
4075 track->SetSigmaY(layer,erry);
4076 track->SetSigmaZ(layer, errz);
4077 track->SetSigmaYZ(layer,covyz);
4078 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4079 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4083 //------------------------------------------------------------------------
4084 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4089 Int_t layer = (index & 0xf0000000) >> 28;
4090 track->SetClIndex(layer, index);
4091 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4092 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4093 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4094 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4098 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
4101 // Take into account the mis-alignment (bring track to cluster plane)
4102 Double_t xTrOrig=track->GetX();
4103 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4104 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4105 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4106 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
4108 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4111 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4112 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4113 c.SetSigmaYZ(track->GetSigmaYZ(layer));
4116 Int_t updated = track->UpdateMI(&c,chi2,index);
4117 // Bring the track back to detector plane in ideal geometry
4118 if (!track->Propagate(xTrOrig)) return 0;
4120 if(!updated) AliDebug(2,"update failed");
4124 //------------------------------------------------------------------------
4125 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4128 //DCA sigmas parameterization
4129 //to be paramterized using external parameters in future
4132 Double_t curv=track->GetC();
4133 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4134 sigmaz = 0.0110+4.37*TMath::Abs(curv);
4136 //------------------------------------------------------------------------
4137 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4140 // Clusters from delta electrons?
4142 Int_t entries = clusterArray->GetEntriesFast();
4143 if (entries<4) return;
4144 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4145 Int_t layer = cluster->GetLayer();
4146 if (layer>1) return;
4148 Int_t ncandidates=0;
4149 Float_t r = (layer>0)? 7:4;
4151 for (Int_t i=0;i<entries;i++){
4152 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4153 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4154 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4155 index[ncandidates] = i; //candidate to belong to delta electron track
4157 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4158 cl0->SetDeltaProbability(1);
4164 for (Int_t i=0;i<ncandidates;i++){
4165 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4166 if (cl0->GetDeltaProbability()>0.8) continue;
4169 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4170 sumy=sumz=sumy2=sumyz=sumw=0.0;
4171 for (Int_t j=0;j<ncandidates;j++){
4173 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4175 Float_t dz = cl0->GetZ()-cl1->GetZ();
4176 Float_t dy = cl0->GetY()-cl1->GetY();
4177 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4178 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4179 y[ncl] = cl1->GetY();
4180 z[ncl] = cl1->GetZ();
4181 sumy+= y[ncl]*weight;
4182 sumz+= z[ncl]*weight;
4183 sumy2+=y[ncl]*y[ncl]*weight;
4184 sumyz+=y[ncl]*z[ncl]*weight;
4189 if (ncl<4) continue;
4190 Float_t det = sumw*sumy2 - sumy*sumy;
4192 if (TMath::Abs(det)>0.01){
4193 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
4194 Float_t k = (sumyz*sumw - sumy*sumz)/det;
4195 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4198 Float_t z0 = sumyz/sumy;
4199 delta = TMath::Abs(cl0->GetZ()-z0);
4202 cl0->SetDeltaProbability(1-20.*delta);
4206 //------------------------------------------------------------------------
4207 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4212 track->UpdateESDtrack(flags);
4213 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4214 if (oldtrack) delete oldtrack;
4215 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4216 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4217 // printf("Problem\n");
4220 //------------------------------------------------------------------------
4221 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4223 // Get nearest upper layer close to the point xr.
4224 // rough approximation
4226 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4227 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4229 for (Int_t i=0;i<6;i++){
4230 if (radius<kRadiuses[i]){
4237 //------------------------------------------------------------------------
4238 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4239 //--------------------------------------------------------------------
4240 // Fill a look-up table with mean material
4241 //--------------------------------------------------------------------
4245 Double_t point1[3],point2[3];
4246 Double_t phi,cosphi,sinphi,z;
4247 // 0-5 layers, 6 pipe, 7-8 shields
4248 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4249 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4251 Int_t ifirst=0,ilast=0;
4252 if(material.Contains("Pipe")) {
4254 } else if(material.Contains("Shields")) {
4256 } else if(material.Contains("Layers")) {
4259 Error("BuildMaterialLUT","Wrong layer name\n");
4262 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4263 Double_t param[5]={0.,0.,0.,0.,0.};
4264 for (Int_t i=0; i<n; i++) {
4265 phi = 2.*TMath::Pi()*gRandom->Rndm();
4266 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4267 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4268 point1[0] = rmin[imat]*cosphi;
4269 point1[1] = rmin[imat]*sinphi;
4271 point2[0] = rmax[imat]*cosphi;
4272 point2[1] = rmax[imat]*sinphi;
4274 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4275 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4277 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4279 fxOverX0Layer[imat] = param[1];
4280 fxTimesRhoLayer[imat] = param[0]*param[4];
4281 } else if(imat==6) {
4282 fxOverX0Pipe = param[1];
4283 fxTimesRhoPipe = param[0]*param[4];
4284 } else if(imat==7) {
4285 fxOverX0Shield[0] = param[1];
4286 fxTimesRhoShield[0] = param[0]*param[4];
4287 } else if(imat==8) {
4288 fxOverX0Shield[1] = param[1];
4289 fxTimesRhoShield[1] = param[0]*param[4];
4293 printf("%s\n",material.Data());
4294 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4295 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4296 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4297 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4298 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4299 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4300 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4301 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4302 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4306 //------------------------------------------------------------------------
4307 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4308 TString direction) {
4309 //-------------------------------------------------------------------
4310 // Propagate beyond beam pipe and correct for material
4311 // (material budget in different ways according to fUseTGeo value)
4312 // Add time if going outward (PropagateTo or PropagateToTGeo)
4313 //-------------------------------------------------------------------
4315 // Define budget mode:
4316 // 0: material from AliITSRecoParam (hard coded)
4317 // 1: material from TGeo in one step (on the fly)
4318 // 2: material from lut
4319 // 3: material from TGeo in one step (same for all hypotheses)
4332 if(fTrackingPhase.Contains("Clusters2Tracks"))
4333 { mode=3; } else { mode=1; }
4336 if(fTrackingPhase.Contains("Clusters2Tracks"))
4337 { mode=3; } else { mode=2; }
4343 if(fTrackingPhase.Contains("Default")) mode=0;
4345 Int_t index=fCurrentEsdTrack;
4347 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4348 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4350 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4352 Double_t xOverX0,x0,lengthTimesMeanDensity;
4356 xOverX0 = AliITSRecoParam::GetdPipe();
4357 x0 = AliITSRecoParam::GetX0Be();
4358 lengthTimesMeanDensity = xOverX0*x0;
4359 lengthTimesMeanDensity *= dir;
4360 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4363 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4366 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4367 xOverX0 = fxOverX0Pipe;
4368 lengthTimesMeanDensity = fxTimesRhoPipe;
4369 lengthTimesMeanDensity *= dir;
4370 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4373 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4374 if(fxOverX0PipeTrks[index]<0) {
4375 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4376 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4377 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4378 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4379 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4382 xOverX0 = fxOverX0PipeTrks[index];
4383 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4384 lengthTimesMeanDensity *= dir;
4385 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4391 //------------------------------------------------------------------------
4392 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4394 TString direction) {
4395 //-------------------------------------------------------------------
4396 // Propagate beyond SPD or SDD shield and correct for material
4397 // (material budget in different ways according to fUseTGeo value)
4398 // Add time if going outward (PropagateTo or PropagateToTGeo)
4399 //-------------------------------------------------------------------
4401 // Define budget mode:
4402 // 0: material from AliITSRecoParam (hard coded)
4403 // 1: material from TGeo in steps of X cm (on the fly)
4404 // X = AliITSRecoParam::GetStepSizeTGeo()
4405 // 2: material from lut
4406 // 3: material from TGeo in one step (same for all hypotheses)
4419 if(fTrackingPhase.Contains("Clusters2Tracks"))
4420 { mode=3; } else { mode=1; }
4423 if(fTrackingPhase.Contains("Clusters2Tracks"))
4424 { mode=3; } else { mode=2; }
4430 if(fTrackingPhase.Contains("Default")) mode=0;
4432 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4434 Int_t shieldindex=0;
4435 if (shield.Contains("SDD")) { // SDDouter
4436 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4438 } else if (shield.Contains("SPD")) { // SPDouter
4439 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4442 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4446 // do nothing if we are already beyond the shield
4447 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4448 if(dir<0 && rTrack > rToGo) return 1; // going outward
4449 if(dir>0 && rTrack < rToGo) return 1; // going inward
4453 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4455 Int_t index=2*fCurrentEsdTrack+shieldindex;
4457 Double_t xOverX0,x0,lengthTimesMeanDensity;
4462 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4463 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4464 lengthTimesMeanDensity = xOverX0*x0;
4465 lengthTimesMeanDensity *= dir;
4466 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4469 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4470 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4473 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4474 xOverX0 = fxOverX0Shield[shieldindex];
4475 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4476 lengthTimesMeanDensity *= dir;
4477 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4480 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4481 if(fxOverX0ShieldTrks[index]<0) {
4482 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4483 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4484 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4485 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4486 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4489 xOverX0 = fxOverX0ShieldTrks[index];
4490 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4491 lengthTimesMeanDensity *= dir;
4492 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4498 //------------------------------------------------------------------------
4499 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4501 Double_t oldGlobXYZ[3],
4502 TString direction) {
4503 //-------------------------------------------------------------------
4504 // Propagate beyond layer and correct for material
4505 // (material budget in different ways according to fUseTGeo value)
4506 // Add time if going outward (PropagateTo or PropagateToTGeo)
4507 //-------------------------------------------------------------------
4509 // Define budget mode:
4510 // 0: material from AliITSRecoParam (hard coded)
4511 // 1: material from TGeo in stepsof X cm (on the fly)
4512 // X = AliITSRecoParam::GetStepSizeTGeo()
4513 // 2: material from lut
4514 // 3: material from TGeo in one step (same for all hypotheses)
4527 if(fTrackingPhase.Contains("Clusters2Tracks"))
4528 { mode=3; } else { mode=1; }
4531 if(fTrackingPhase.Contains("Clusters2Tracks"))
4532 { mode=3; } else { mode=2; }
4538 if(fTrackingPhase.Contains("Default")) mode=0;
4540 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4542 Double_t r=fgLayers[layerindex].GetR();
4543 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4545 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4547 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4549 Int_t index=6*fCurrentEsdTrack+layerindex;
4552 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4555 // back before material (no correction)
4557 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4558 if (!t->GetLocalXat(rOld,xOld)) return 0;
4559 if (!t->Propagate(xOld)) return 0;
4563 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4564 lengthTimesMeanDensity = xOverX0*x0;
4565 lengthTimesMeanDensity *= dir;
4566 // Bring the track beyond the material
4567 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4570 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4571 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4574 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4575 xOverX0 = fxOverX0Layer[layerindex];
4576 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4577 lengthTimesMeanDensity *= dir;
4578 // Bring the track beyond the material
4579 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4582 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4583 if(fxOverX0LayerTrks[index]<0) {
4584 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4585 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4586 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4587 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4588 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4589 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4592 xOverX0 = fxOverX0LayerTrks[index];
4593 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4594 lengthTimesMeanDensity *= dir;
4595 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4602 //------------------------------------------------------------------------
4603 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4604 //-----------------------------------------------------------------
4605 // Initialize LUT for storing material for each prolonged track
4606 //-----------------------------------------------------------------
4607 fxOverX0PipeTrks = new Float_t[ntracks];
4608 fxTimesRhoPipeTrks = new Float_t[ntracks];
4609 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4610 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4611 fxOverX0LayerTrks = new Float_t[ntracks*6];
4612 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4614 for(Int_t i=0; i<ntracks; i++) {
4615 fxOverX0PipeTrks[i] = -1.;
4616 fxTimesRhoPipeTrks[i] = -1.;
4618 for(Int_t j=0; j<ntracks*2; j++) {
4619 fxOverX0ShieldTrks[j] = -1.;
4620 fxTimesRhoShieldTrks[j] = -1.;
4622 for(Int_t k=0; k<ntracks*6; k++) {
4623 fxOverX0LayerTrks[k] = -1.;
4624 fxTimesRhoLayerTrks[k] = -1.;
4631 //------------------------------------------------------------------------
4632 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4633 //-----------------------------------------------------------------
4634 // Delete LUT for storing material for each prolonged track
4635 //-----------------------------------------------------------------
4636 if(fxOverX0PipeTrks) {
4637 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4639 if(fxOverX0ShieldTrks) {
4640 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4643 if(fxOverX0LayerTrks) {
4644 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4646 if(fxTimesRhoPipeTrks) {
4647 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4649 if(fxTimesRhoShieldTrks) {
4650 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4652 if(fxTimesRhoLayerTrks) {
4653 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4657 //------------------------------------------------------------------------
4658 void AliITStrackerMI::SetForceSkippingOfLayer() {
4659 //-----------------------------------------------------------------
4660 // Check if we are forced to skip layers
4661 // either we set to skip them in RecoParam
4662 // or they were off during data-taking
4663 //-----------------------------------------------------------------
4665 const AliEventInfo *eventInfo = GetEventInfo();
4667 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4668 fForceSkippingOfLayer[l] = 0;
4670 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4674 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4675 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4677 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4678 } else if(l==2 || l==3) {
4679 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4681 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4687 //------------------------------------------------------------------------
4688 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4689 Int_t ilayer,Int_t idet) const {
4690 //-----------------------------------------------------------------
4691 // This method is used to decide whether to allow a prolongation
4692 // without clusters, because we want to skip the layer.
4693 // In this case the return value is > 0:
4694 // return 1: the user requested to skip a layer
4695 // return 2: track outside z acceptance
4696 //-----------------------------------------------------------------
4698 if (ForceSkippingOfLayer(ilayer)) return 1;
4700 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4702 if (idet<0 && // out in z
4703 ilayer>innerLayCanSkip &&
4704 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4705 // check if track will cross SPD outer layer
4706 Double_t phiAtSPD2,zAtSPD2;
4707 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4708 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4710 return 2; // always allow skipping, changed on 05.11.2009
4715 //------------------------------------------------------------------------
4716 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4717 Int_t ilayer,Int_t idet,
4718 Double_t dz,Double_t dy,
4719 Bool_t noClusters) const {
4720 //-----------------------------------------------------------------
4721 // This method is used to decide whether to allow a prolongation
4722 // without clusters, because there is a dead zone in the road.
4723 // In this case the return value is > 0:
4724 // return 1: dead zone at z=0,+-7cm in SPD
4725 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4726 // return 2: all road is "bad" (dead or noisy) from the OCDB
4727 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4728 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4729 //-----------------------------------------------------------------
4731 // check dead zones at z=0,+-7cm in the SPD
4732 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4733 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4734 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4735 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4736 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4737 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4738 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4739 for (Int_t i=0; i<3; i++)
4740 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4741 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4742 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4746 // check bad zones from OCDB
4747 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4749 if (idet<0) return 0;
4751 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4754 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4755 if (ilayer==0 || ilayer==1) { // ---------- SPD
4757 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4759 detSizeFactorX *= 2.;
4760 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4763 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4764 if (detType==2) segm->SetLayer(ilayer+1);
4765 Float_t detSizeX = detSizeFactorX*segm->Dx();
4766 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4768 // check if the road overlaps with bad chips
4770 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4771 Float_t zlocmin = zloc-dz;
4772 Float_t zlocmax = zloc+dz;
4773 Float_t xlocmin = xloc-dy;
4774 Float_t xlocmax = xloc+dy;
4775 Int_t chipsInRoad[100];
4777 // check if road goes out of detector
4778 Bool_t touchNeighbourDet=kFALSE;
4779 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4780 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4781 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4782 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4783 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));
4785 // check if this detector is bad
4787 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4788 if(!touchNeighbourDet) {
4789 return 2; // all detectors in road are bad
4791 return 3; // at least one is bad
4795 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4796 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4797 if (!nChipsInRoad) return 0;
4799 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4800 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4801 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4802 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4803 if (det.IsChipBad(chipsInRoad[iCh])) {
4811 if(!touchNeighbourDet) {
4812 AliDebug(2,"all bad in road");
4813 return 2; // all chips in road are bad
4815 return 3; // at least a bad chip in road
4820 AliDebug(2,"at least a bad in road");
4821 return 3; // at least a bad chip in road
4825 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4826 || !noClusters) return 0;
4828 // There are no clusters in road: check if there is at least
4829 // a bad SPD pixel or SDD anode or SSD strips on both sides
4831 Int_t idetInITS=idet;
4832 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4834 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4835 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4838 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4842 //------------------------------------------------------------------------
4843 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4844 const AliITStrackMI *track,
4845 Float_t &xloc,Float_t &zloc) const {
4846 //-----------------------------------------------------------------
4847 // Gives position of track in local module ref. frame
4848 //-----------------------------------------------------------------
4853 if(idet<0) return kTRUE; // track out of z acceptance of layer
4855 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4857 Int_t lad = Int_t(idet/ndet) + 1;
4859 Int_t det = idet - (lad-1)*ndet + 1;
4861 Double_t xyzGlob[3],xyzLoc[3];
4863 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4864 // take into account the misalignment: xyz at real detector plane
4865 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4867 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4869 xloc = (Float_t)xyzLoc[0];
4870 zloc = (Float_t)xyzLoc[2];
4874 //------------------------------------------------------------------------
4875 //------------------------------------------------------------------------
4876 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4878 // Method to be optimized further:
4879 // Aim: decide whether a track can be used for PlaneEff evaluation
4880 // the decision is taken based on the track quality at the layer under study
4881 // no information on the clusters on this layer has to be used
4882 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4883 // the cut is done on number of sigmas from the boundaries
4885 // Input: Actual track, layer [0,5] under study
4887 // Return: kTRUE if this is a good track
4889 // it will apply a pre-selection to obtain good quality tracks.
4890 // Here also you will have the possibility to put a control on the
4891 // impact point of the track on the basic block, in order to exclude border regions
4892 // this will be done by calling a proper method of the AliITSPlaneEff class.
4894 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4895 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4897 Int_t index[AliITSgeomTGeo::kNLayers];
4899 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4901 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4902 index[k]=clusters[k];
4906 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4907 AliITSlayer &layer=fgLayers[ilayer];
4908 Double_t r=layer.GetR();
4909 AliITStrackMI tmp(*track);
4911 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4912 Int_t ncl_out=0; Int_t ncl_in=0;
4913 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4914 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4915 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4916 // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4917 if(index[lay]>=0)ncl_out++;
4919 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4920 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4921 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4922 if (index[lay]>=0) ncl_in++;
4924 Int_t ncl=ncl_out+ncl_in;
4925 Bool_t nextout = kFALSE;
4926 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4927 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4928 Bool_t nextin = kFALSE;
4929 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4930 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4931 // maximum number of missing clusters allowed in outermost layers
4932 if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4934 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4935 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4937 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4938 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4939 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4940 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4944 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4945 Int_t idet=layer.FindDetectorIndex(phi,z);
4946 if(idet<0) { AliInfo(Form("cannot find detector"));
4949 // here check if it has good Chi Square.
4951 //propagate to the intersection with the detector plane
4952 const AliITSdetector &det=layer.GetDetector(idet);
4953 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4957 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4958 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4959 if(key>fPlaneEff->Nblock()) return kFALSE;
4960 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4961 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4963 // DEFINITION OF SEARCH ROAD FOR accepting a track
4965 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4966 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4967 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4968 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4969 // done for RecPoints
4971 // exclude tracks at boundary between detectors
4972 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4973 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4974 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4975 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4976 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4977 if ( (locx-dx < blockXmn+boundaryWidth) ||
4978 (locx+dx > blockXmx-boundaryWidth) ||
4979 (locz-dz < blockZmn+boundaryWidth) ||
4980 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4983 //------------------------------------------------------------------------
4984 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4986 // This Method has to be optimized! For the time-being it uses the same criteria
4987 // as those used in the search of extra clusters for overlapping modules.
4989 // Method Purpose: estabilish whether a track has produced a recpoint or not
4990 // in the layer under study (For Plane efficiency)
4992 // inputs: AliITStrackMI* track (pointer to a usable track)
4994 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4995 // traversed by this very track. In details:
4996 // - if a cluster can be associated to the track then call
4997 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4998 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5001 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5002 AliITSlayer &layer=fgLayers[ilayer];
5003 Double_t r=layer.GetR();
5004 AliITStrackMI tmp(*track);
5008 if (!tmp.GetPhiZat(r,phi,z)) return;
5009 Int_t idet=layer.FindDetectorIndex(phi,z);
5011 if(idet<0) { AliInfo(Form("cannot find detector"));
5015 //propagate to the intersection with the detector plane
5016 const AliITSdetector &det=layer.GetDetector(idet);
5017 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5021 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5023 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5024 TMath::Sqrt(tmp.GetSigmaZ2() +
5025 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5026 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5027 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5028 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5029 TMath::Sqrt(tmp.GetSigmaY2() +
5030 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5031 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5032 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5034 // road in global (rphi,z) [i.e. in tracking ref. system]
5035 Double_t zmin = tmp.GetZ() - dz;
5036 Double_t zmax = tmp.GetZ() + dz;
5037 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5038 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5040 // select clusters in road
5041 layer.SelectClusters(zmin,zmax,ymin,ymax);
5043 // Define criteria for track-cluster association
5044 Double_t msz = tmp.GetSigmaZ2() +
5045 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5046 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5047 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5048 Double_t msy = tmp.GetSigmaY2() +
5049 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5050 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5051 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5052 if (tmp.GetConstrain()) {
5053 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5054 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5056 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5057 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5059 msz = 1./msz; // 1/RoadZ^2
5060 msy = 1./msy; // 1/RoadY^2
5063 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5065 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5066 //Double_t tolerance=0.2;
5067 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5068 idetc = cl->GetDetectorIndex();
5069 if(idet!=idetc) continue;
5070 //Int_t ilay = cl->GetLayer();
5072 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5073 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5075 Double_t chi2=tmp.GetPredictedChi2(cl);
5076 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5080 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5082 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5083 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5084 if(key>fPlaneEff->Nblock()) return;
5085 Bool_t found=kFALSE;
5088 while ((cl=layer.GetNextCluster(clidx))!=0) {
5089 idetc = cl->GetDetectorIndex();
5090 if(idet!=idetc) continue;
5091 // here real control to see whether the cluster can be associated to the track.
5092 // cluster not associated to track
5093 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5094 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5095 // calculate track-clusters chi2
5096 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5097 // in particular, the error associated to the cluster
5098 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5100 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5102 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5103 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5104 // track->SetExtraModule(ilayer,idetExtra);
5106 if(!fPlaneEff->UpDatePlaneEff(found,key))
5107 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5109 // this for FO efficiency studies (only for SPD) //
5110 UInt_t keyFO=999999;
5111 Bool_t foundFO=kFALSE;
5112 if(ilayer<2){ //ONLY SPD layers for FastOr studies
5113 TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5114 Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5115 if(!fSPDChipIntPlaneEff[key]){
5116 AliITSPlaneEffSPD spd;
5117 keyFO = spd.SwitchChipKeyNumbering(key);
5118 if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5119 keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5120 if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5121 AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5124 if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5125 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5131 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5132 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5133 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5134 Int_t cltype[2]={-999,-999};
5137 Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5141 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5142 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5145 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5146 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5147 cltype[0]=layer.GetCluster(ci)->GetNy();
5148 cltype[1]=layer.GetCluster(ci)->GetNz();
5150 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5151 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5152 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5153 // It is computed properly by calling the method
5154 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5156 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5157 //if (tmp.PropagateTo(x,0.,0.)) {
5158 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5159 AliCluster c(*layer.GetCluster(ci));
5160 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5161 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5162 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5163 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5164 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5167 // Compute the angles between the track and the module
5168 // compute the angle "in phi direction", i.e. the angle in the transverse plane
5169 // between the normal to the module and the projection (in the transverse plane) of the
5171 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5172 Float_t tgl = tmp.GetTgl();
5173 Float_t phitr = tmp.GetSnp();
5174 phitr = TMath::ASin(phitr);
5175 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5177 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5178 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5180 alpha = tmp.GetAlpha();
5181 Double_t phiglob = alpha+phitr;
5183 p[0] = TMath::Cos(phiglob);
5184 p[1] = TMath::Sin(phiglob);
5186 TVector3 pvec(p[0],p[1],p[2]);
5187 TVector3 normvec(rot[1],rot[4],rot[7]);
5188 Double_t angle = pvec.Angle(normvec);
5190 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5191 angle *= 180./TMath::Pi();
5194 TVector3 pt(p[0],p[1],0);
5195 TVector3 normt(rot[1],rot[4],0);
5196 Double_t anglet = pt.Angle(normt);
5198 Double_t phiPt = TMath::ATan2(p[1],p[0]);
5199 if(phiPt<0)phiPt+=2.*TMath::Pi();
5200 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5201 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5202 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5203 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
5204 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5205 anglet *= 180./TMath::Pi();
5207 AngleModTrack[2]=(Float_t) angle;
5208 AngleModTrack[0]=(Float_t) anglet;
5209 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5210 AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5211 AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5212 AngleModTrack[1]*=180./TMath::Pi(); // in degree
5214 fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
5216 // For FO efficiency studies of SPD
5217 if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,AngleModTrack);
5219 if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5223 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5225 // find the MC cluster for the label
5226 return fgLayers[lr].FindClusterForLabel(label,store);
5230 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5232 // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5233 strncpy(patt,"......",6);
5235 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5238 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5239 Int_t cindex = track->GetClusterIndex(i);
5240 Int_t l=(cindex & 0xf0000000) >> 28;
5241 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5243 for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5244 patt[l] = isWrong ? 'f':'c';
5250 //------------------------------------------------------------------------
5251 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5253 //--------------------------------------------------------------------
5256 for (int ic=0;ic<fN;ic++) {
5257 const AliITSRecPoint *cl = GetCluster(ic);
5258 for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5260 if (store) store[nfound] = ic;