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(*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)) continue;
738 fTrackToFollow.SetLabel(t.GetLabel());
739 //fTrackToFollow.CookdEdx();
740 CookLabel(&fTrackToFollow,0.); //For comparison only
741 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
742 //UseClusters(&fTrackToFollow);
747 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
749 fTrackingPhase="Default";
753 //------------------------------------------------------------------------
754 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
755 //--------------------------------------------------------------------
756 // This functions refits ITS tracks using the
757 // "inward propagated" TPC tracks
758 // The clusters must be loaded !
759 //--------------------------------------------------------------------
760 fTrackingPhase="RefitInward";
762 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
764 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
765 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
767 Int_t nentr=event->GetNumberOfTracks();
768 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
770 // only for PlaneEff and in case of SPD (for FO studies)
771 if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
772 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 &&
773 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) {
774 for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
778 for (Int_t i=0; i<nentr; i++) {
779 AliESDtrack *esd=event->GetTrack(i);
781 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
782 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
783 if (esd->GetStatus()&AliESDtrack::kTPCout)
784 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
786 AliITStrackMI *t = new AliITStrackMI(*esd);
788 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
789 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
794 ResetTrackToFollow(*t);
795 fTrackToFollow.ResetClusters();
797 // ITS standalone tracks
798 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
799 fTrackToFollow.ResetCovariance(10.);
800 // protection for loopers that can have parameters screwed up
801 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
802 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
809 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
810 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
812 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
813 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
814 AliDebug(2," refit OK");
815 fTrackToFollow.SetLabel(t->GetLabel());
816 // fTrackToFollow.CookdEdx();
817 CookdEdx(&fTrackToFollow);
819 CookLabel(&fTrackToFollow,0.0); //For comparison only
822 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
823 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
824 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
825 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
826 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
827 Double_t r[3]={0.,0.,0.};
829 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
836 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
838 fTrackingPhase="Default";
842 //------------------------------------------------------------------------
843 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
844 //--------------------------------------------------------------------
845 // Return pointer to a given cluster
846 //--------------------------------------------------------------------
847 Int_t l=(index & 0xf0000000) >> 28;
848 Int_t c=(index & 0x0fffffff) >> 00;
849 return fgLayers[l].GetCluster(c);
851 //------------------------------------------------------------------------
852 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
853 //--------------------------------------------------------------------
854 // Get track space point with index i
855 //--------------------------------------------------------------------
857 Int_t l=(index & 0xf0000000) >> 28;
858 Int_t c=(index & 0x0fffffff) >> 00;
859 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
860 Int_t idet = cl->GetDetectorIndex();
864 cl->GetGlobalXYZ(xyz);
865 cl->GetGlobalCov(cov);
867 p.SetCharge(cl->GetQ());
868 p.SetDriftTime(cl->GetDriftTime());
869 p.SetChargeRatio(cl->GetChargeRatio());
870 p.SetClusterType(cl->GetClusterType());
871 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
874 iLayer = AliGeomManager::kSPD1;
877 iLayer = AliGeomManager::kSPD2;
880 iLayer = AliGeomManager::kSDD1;
883 iLayer = AliGeomManager::kSDD2;
886 iLayer = AliGeomManager::kSSD1;
889 iLayer = AliGeomManager::kSSD2;
892 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
895 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
896 p.SetVolumeID((UShort_t)volid);
899 //------------------------------------------------------------------------
900 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
901 AliTrackPoint& p, const AliESDtrack *t) {
902 //--------------------------------------------------------------------
903 // Get track space point with index i
904 // (assign error estimated during the tracking)
905 //--------------------------------------------------------------------
907 Int_t l=(index & 0xf0000000) >> 28;
908 Int_t c=(index & 0x0fffffff) >> 00;
909 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
910 Int_t idet = cl->GetDetectorIndex();
912 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
914 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
916 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
917 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
918 Double_t alpha = t->GetAlpha();
919 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
920 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
921 phi += alpha-det.GetPhi();
922 Float_t tgphi = TMath::Tan(phi);
924 Float_t tgl = t->GetTgl(); // tgl about const along track
925 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
927 Float_t errtrky,errtrkz,covyz;
928 Bool_t addMisalErr=kFALSE;
929 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
933 cl->GetGlobalXYZ(xyz);
934 // cl->GetGlobalCov(cov);
935 Float_t pos[3] = {0.,0.,0.};
936 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
937 tmpcl.GetGlobalCov(cov);
940 p.SetCharge(cl->GetQ());
941 p.SetDriftTime(cl->GetDriftTime());
942 p.SetChargeRatio(cl->GetChargeRatio());
943 p.SetClusterType(cl->GetClusterType());
945 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
948 iLayer = AliGeomManager::kSPD1;
951 iLayer = AliGeomManager::kSPD2;
954 iLayer = AliGeomManager::kSDD1;
957 iLayer = AliGeomManager::kSDD2;
960 iLayer = AliGeomManager::kSSD1;
963 iLayer = AliGeomManager::kSSD2;
966 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
969 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
971 p.SetVolumeID((UShort_t)volid);
974 //------------------------------------------------------------------------
975 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
977 //--------------------------------------------------------------------
978 // Follow prolongation tree
979 //--------------------------------------------------------------------
981 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
982 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
985 AliESDtrack * esd = otrack->GetESDtrack();
986 if (esd->GetV0Index(0)>0) {
987 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
988 // mapping of ESD track is different as ITS track in Containers
989 // Need something more stable
990 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
991 for (Int_t i=0;i<3;i++){
992 Int_t index = esd->GetV0Index(i);
994 AliESDv0 * vertex = fEsd->GetV0(index);
995 if (vertex->GetStatus()<0) continue; // rejected V0
997 if (esd->GetSign()>0) {
998 vertex->SetIndex(0,esdindex);
1000 vertex->SetIndex(1,esdindex);
1004 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1006 bestarray = new TObjArray(5);
1007 bestarray->SetOwner();
1008 fBestHypothesys.AddAt(bestarray,esdindex);
1012 //setup tree of the prolongations
1014 const int kMaxTr = 100; //RS
1015 static AliITStrackMI tracks[7][kMaxTr];
1016 AliITStrackMI *currenttrack;
1017 static AliITStrackMI currenttrack1;
1018 static AliITStrackMI currenttrack2;
1019 static AliITStrackMI backuptrack;
1021 Int_t nindexes[7][kMaxTr];
1022 Float_t normalizedchi2[kMaxTr];
1023 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1024 otrack->SetNSkipped(0);
1025 new (&(tracks[6][0])) AliITStrackMI(*otrack);
1027 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1028 Int_t modstatus = 1; // found
1032 // follow prolongations
1033 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1034 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1037 AliITSlayer &layer=fgLayers[ilayer];
1038 Double_t r = layer.GetR();
1044 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1045 // printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1047 if (ntracks[ilayer]>=kMaxTr) break;
1048 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1049 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1050 if (ntracks[ilayer]>15+ilayer){
1051 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1052 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1055 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1057 // material between SSD and SDD, SDD and SPD
1059 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1061 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1065 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1066 Int_t idet=layer.FindDetectorIndex(phi,z);
1068 Double_t trackGlobXYZ1[3];
1069 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1071 // Get the budget to the primary vertex for the current track being prolonged
1072 Double_t budgetToPrimVertex = 0;
1073 double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1076 if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1077 else budgetToPrimVertex = GetEffectiveThickness();
1079 // check if we allow a prolongation without point
1080 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1082 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1083 // propagate to the layer radius
1084 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1085 if(!vtrack->Propagate(xToGo)) continue;
1086 // apply correction for material of the current layer
1087 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1088 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1089 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1090 vtrack->SetClIndex(ilayer,-1);
1091 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1092 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1093 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1095 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1097 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1098 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1104 // track outside layer acceptance in z
1105 if (idet<0) continue;
1107 //propagate to the intersection with the detector plane
1108 const AliITSdetector &det=layer.GetDetector(idet);
1109 new(¤ttrack2) AliITStrackMI(currenttrack1);
1110 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1111 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1112 currenttrack1.SetDetectorIndex(idet);
1113 currenttrack2.SetDetectorIndex(idet);
1114 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1117 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1119 // road in global (rphi,z) [i.e. in tracking ref. system]
1120 Double_t zmin,zmax,ymin,ymax;
1121 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1123 // select clusters in road
1124 layer.SelectClusters(zmin,zmax,ymin,ymax);
1125 //********************
1127 // Define criteria for track-cluster association
1128 Double_t msz = currenttrack1.GetSigmaZ2() +
1129 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1130 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1131 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1132 Double_t msy = currenttrack1.GetSigmaY2() +
1133 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1134 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1135 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1137 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1138 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1140 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1141 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1143 msz = 1./msz; // 1/RoadZ^2
1144 msy = 1./msy; // 1/RoadY^2
1148 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1150 const AliITSRecPoint *cl=0;
1152 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1153 Bool_t deadzoneSPD=kFALSE;
1154 currenttrack = ¤ttrack1;
1156 // check if the road contains a dead zone
1157 Bool_t noClusters = kFALSE;
1158 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1159 if (noClusters) AliDebug(2,"no clusters in road");
1160 Double_t dz=0.5*(zmax-zmin);
1161 Double_t dy=0.5*(ymax-ymin);
1162 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1163 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1164 // create a prolongation without clusters (check also if there are no clusters in the road)
1167 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1168 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1169 updatetrack->SetClIndex(ilayer,-1);
1171 modstatus = 5; // no cls in road
1172 } else if (dead==1) {
1173 modstatus = 7; // holes in z in SPD
1174 } else if (dead==2 || dead==3 || dead==4) {
1175 modstatus = 2; // dead from OCDB
1177 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1178 // apply correction for material of the current layer
1179 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1180 if (constrain) { // apply vertex constrain
1181 updatetrack->SetConstrain(constrain);
1182 Bool_t isPrim = kTRUE;
1183 if (ilayer<4) { // check that it's close to the vertex
1184 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1185 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1186 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1187 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1188 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1190 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1192 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1193 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1196 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1198 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1199 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1201 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1202 updatetrack->SetDeadZoneProbability(ilayer,1.);
1203 } else if (dead==4) { // at least a single dead channel from OCDB
1204 updatetrack->SetDeadZoneProbability(ilayer,0.);
1211 // loop over clusters in the road
1212 while ((cl=layer.GetNextCluster(clidx))!=0) {
1213 if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters
1214 Bool_t changedet =kFALSE;
1215 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1216 Int_t idetc=cl->GetDetectorIndex();
1218 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1219 // take into account misalignment (bring track to real detector plane)
1220 Double_t xTrOrig = currenttrack->GetX();
1221 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1222 // a first cut on track-cluster distance
1223 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1224 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1225 { // cluster not associated to track
1226 AliDebug(2,"not associated");
1227 // MvL: added here as well
1228 // bring track back to ideal detector plane
1229 currenttrack->Propagate(xTrOrig);
1232 // bring track back to ideal detector plane
1233 if (!currenttrack->Propagate(xTrOrig)) continue;
1234 } else { // have to move track to cluster's detector
1235 const AliITSdetector &detc=layer.GetDetector(idetc);
1236 // a first cut on track-cluster distance
1238 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1239 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1240 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1241 continue; // cluster not associated to track
1243 new (&backuptrack) AliITStrackMI(currenttrack2);
1245 currenttrack =¤ttrack2;
1246 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1247 new (currenttrack) AliITStrackMI(backuptrack);
1251 currenttrack->SetDetectorIndex(idetc);
1252 // Get again the budget to the primary vertex
1253 // for the current track being prolonged, if had to change detector
1254 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1257 // calculate track-clusters chi2
1258 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1260 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1261 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1262 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1263 if (ntracks[ilayer]>=kMaxTr) continue;
1264 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1265 updatetrack->SetClIndex(ilayer,-1);
1266 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1268 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1269 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1270 AliDebug(2,"update failed");
1273 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1274 modstatus = 1; // found
1275 } else { // virtual cluster in dead zone
1276 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1277 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1278 modstatus = 7; // holes in z in SPD
1282 Float_t xlocnewdet,zlocnewdet;
1283 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1284 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1287 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1289 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1291 // apply correction for material of the current layer
1292 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1294 if (constrain) { // apply vertex constrain
1295 updatetrack->SetConstrain(constrain);
1296 Bool_t isPrim = kTRUE;
1297 if (ilayer<4) { // check that it's close to the vertex
1298 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1299 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1300 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1301 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1302 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1304 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1306 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1307 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1309 } //apply vertex constrain
1311 } // create new hypothesis
1313 AliDebug(2,"chi2 too large");
1316 } // loop over possible prolongations
1318 // allow one prolongation without clusters
1319 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1320 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1321 // apply correction for material of the current layer
1322 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1323 vtrack->SetClIndex(ilayer,-1);
1324 modstatus = 3; // skipped
1325 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1326 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1328 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1329 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1331 vtrack->IncrementNSkipped();
1336 } // loop over tracks in layer ilayer+1
1338 //loop over track candidates for the current layer
1344 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1345 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1346 if (normalizedchi2[itrack] <
1347 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1351 if (constrain) { // constrain
1352 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1354 } else { // !constrain
1355 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1360 // sort tracks by increasing normalized chi2
1361 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1362 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1363 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1364 // if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1365 if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9);
1366 } // end loop over layers
1370 // Now select tracks to be kept
1372 Int_t max = constrain ? 20 : 5;
1374 // tracks that reach layer 0 (SPD inner)
1375 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1376 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1377 if (track.GetNumberOfClusters()<2) continue;
1378 if (!constrain && track.GetNormChi2(0) >
1379 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1382 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1385 // tracks that reach layer 1 (SPD outer)
1386 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1387 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1388 if (track.GetNumberOfClusters()<4) continue;
1389 if (!constrain && track.GetNormChi2(1) >
1390 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1391 if (constrain) track.IncrementNSkipped();
1393 track.SetD(0,track.GetD(GetX(),GetY()));
1394 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1395 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1396 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1399 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1402 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1404 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1405 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1406 if (track.GetNumberOfClusters()<3) continue;
1407 if (track.GetNormChi2(2) >
1408 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1409 track.SetD(0,track.GetD(GetX(),GetY()));
1410 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1411 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1412 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1414 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1420 // register best track of each layer - important for V0 finder
1422 for (Int_t ilayer=0;ilayer<5;ilayer++){
1423 if (ntracks[ilayer]==0) continue;
1424 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1425 if (track.GetNumberOfClusters()<1) continue;
1426 CookLabel(&track,0);
1427 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1431 // update TPC V0 information
1433 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1434 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1435 for (Int_t i=0;i<3;i++){
1436 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1437 if (index==0) break;
1438 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1439 if (vertex->GetStatus()<0) continue; // rejected V0
1441 if (otrack->GetSign()>0) {
1442 vertex->SetIndex(0,esdindex);
1445 vertex->SetIndex(1,esdindex);
1447 //find nearest layer with track info
1448 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1449 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1450 Int_t nearest = nearestold;
1451 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1452 if (ntracks[nearest]==0){
1457 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1458 if (nearestold<5&&nearest<5){
1459 Bool_t accept = track.GetNormChi2(nearest)<10;
1461 if (track.GetSign()>0) {
1462 vertex->SetParamP(track);
1463 vertex->Update(fprimvertex);
1464 //vertex->SetIndex(0,track.fESDtrack->GetID());
1465 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1467 vertex->SetParamN(track);
1468 vertex->Update(fprimvertex);
1469 //vertex->SetIndex(1,track.fESDtrack->GetID());
1470 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1472 vertex->SetStatus(vertex->GetStatus()+1);
1474 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1481 //------------------------------------------------------------------------
1482 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1484 //--------------------------------------------------------------------
1487 return fgLayers[layer];
1489 //------------------------------------------------------------------------
1490 AliITStrackerMI::AliITSlayer::AliITSlayer():
1520 //--------------------------------------------------------------------
1521 //default AliITSlayer constructor
1522 //--------------------------------------------------------------------
1523 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1524 fClusterWeight[i]=0;
1525 fClusterTracks[0][i]=-1;
1526 fClusterTracks[1][i]=-1;
1527 fClusterTracks[2][i]=-1;
1528 fClusterTracks[3][i]=-1;
1535 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1536 for (Int_t j1=0; j1<6; j1++) {
1537 fClusters5[j1][j]=0;
1538 fClusterIndex5[j1][j]=-1;
1547 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1548 for (Int_t j1=0; j1<11; j1++) {
1549 fClusters10[j1][j]=0;
1550 fClusterIndex10[j1][j]=-1;
1559 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1560 for (Int_t j1=0; j1<21; j1++) {
1561 fClusters20[j1][j]=0;
1562 fClusterIndex20[j1][j]=-1;
1570 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1575 //------------------------------------------------------------------------
1576 AliITStrackerMI::AliITSlayer::
1577 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1606 //--------------------------------------------------------------------
1607 //main AliITSlayer constructor
1608 //--------------------------------------------------------------------
1609 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1610 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1612 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1613 fClusterWeight[i]=0;
1614 fClusterTracks[0][i]=-1;
1615 fClusterTracks[1][i]=-1;
1616 fClusterTracks[2][i]=-1;
1617 fClusterTracks[3][i]=-1;
1625 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1626 for (Int_t j1=0; j1<6; j1++) {
1627 fClusters5[j1][j]=0;
1628 fClusterIndex5[j1][j]=-1;
1637 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1638 for (Int_t j1=0; j1<11; j1++) {
1639 fClusters10[j1][j]=0;
1640 fClusterIndex10[j1][j]=-1;
1649 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1650 for (Int_t j1=0; j1<21; j1++) {
1651 fClusters20[j1][j]=0;
1652 fClusterIndex20[j1][j]=-1;
1660 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1666 //------------------------------------------------------------------------
1667 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1669 fPhiOffset(layer.fPhiOffset),
1670 fNladders(layer.fNladders),
1671 fZOffset(layer.fZOffset),
1672 fNdetectors(layer.fNdetectors),
1673 fDetectors(layer.fDetectors),
1678 fClustersCs(layer.fClustersCs),
1679 fClusterIndexCs(layer.fClusterIndexCs),
1683 fCurrentSlice(layer.fCurrentSlice),
1691 fAccepted(layer.fAccepted),
1693 fMaxSigmaClY(layer.fMaxSigmaClY),
1694 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1695 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1700 //------------------------------------------------------------------------
1701 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1702 //--------------------------------------------------------------------
1703 // AliITSlayer destructor
1704 //--------------------------------------------------------------------
1705 delete [] fDetectors;
1706 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1707 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1708 fClusterWeight[i]=0;
1709 fClusterTracks[0][i]=-1;
1710 fClusterTracks[1][i]=-1;
1711 fClusterTracks[2][i]=-1;
1712 fClusterTracks[3][i]=-1;
1715 //------------------------------------------------------------------------
1716 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1717 //--------------------------------------------------------------------
1718 // This function removes loaded clusters
1719 //--------------------------------------------------------------------
1720 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1721 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1722 fClusterWeight[i]=0;
1723 fClusterTracks[0][i]=-1;
1724 fClusterTracks[1][i]=-1;
1725 fClusterTracks[2][i]=-1;
1726 fClusterTracks[3][i]=-1;
1732 //------------------------------------------------------------------------
1733 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1734 //--------------------------------------------------------------------
1735 // This function reset weights of the clusters
1736 //--------------------------------------------------------------------
1737 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1738 fClusterWeight[i]=0;
1739 fClusterTracks[0][i]=-1;
1740 fClusterTracks[1][i]=-1;
1741 fClusterTracks[2][i]=-1;
1742 fClusterTracks[3][i]=-1;
1744 for (Int_t i=0; i<fN;i++) {
1745 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1746 if (cl&&cl->IsUsed()) cl->Use();
1750 //------------------------------------------------------------------------
1751 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1752 //--------------------------------------------------------------------
1753 // This function calculates the road defined by the cluster density
1754 //--------------------------------------------------------------------
1756 for (Int_t i=0; i<fN; i++) {
1757 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1759 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1761 //------------------------------------------------------------------------
1762 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1763 //--------------------------------------------------------------------
1764 //This function adds a cluster to this layer
1765 //--------------------------------------------------------------------
1766 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1772 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1774 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1775 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1776 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1777 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1778 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1779 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1782 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1783 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1784 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1785 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1789 //------------------------------------------------------------------------
1790 void AliITStrackerMI::AliITSlayer::SortClusters()
1795 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1796 Float_t *z = new Float_t[fN];
1797 Int_t * index = new Int_t[fN];
1799 fMaxSigmaClY=0.; //AD
1800 fMaxSigmaClZ=0.; //AD
1802 for (Int_t i=0;i<fN;i++){
1803 z[i] = fClusters[i]->GetZ();
1804 // save largest errors in y and z for this layer
1805 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1806 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1808 TMath::Sort(fN,z,index,kFALSE);
1809 for (Int_t i=0;i<fN;i++){
1810 clusters[i] = fClusters[index[i]];
1813 for (Int_t i=0;i<fN;i++){
1814 fClusters[i] = clusters[i];
1815 fZ[i] = fClusters[i]->GetZ();
1816 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1817 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1818 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1828 for (Int_t i=0;i<fN;i++){
1829 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1830 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1831 fClusterIndex[i] = i;
1835 fDy5 = (fYB[1]-fYB[0])/5.;
1836 fDy10 = (fYB[1]-fYB[0])/10.;
1837 fDy20 = (fYB[1]-fYB[0])/20.;
1838 for (Int_t i=0;i<6;i++) fN5[i] =0;
1839 for (Int_t i=0;i<11;i++) fN10[i]=0;
1840 for (Int_t i=0;i<21;i++) fN20[i]=0;
1842 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;}
1843 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;}
1844 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;}
1847 for (Int_t i=0;i<fN;i++)
1848 for (Int_t irot=-1;irot<=1;irot++){
1849 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1851 for (Int_t slice=0; slice<6;slice++){
1852 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1853 fClusters5[slice][fN5[slice]] = fClusters[i];
1854 fY5[slice][fN5[slice]] = curY;
1855 fZ5[slice][fN5[slice]] = fZ[i];
1856 fClusterIndex5[slice][fN5[slice]]=i;
1861 for (Int_t slice=0; slice<11;slice++){
1862 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1863 fClusters10[slice][fN10[slice]] = fClusters[i];
1864 fY10[slice][fN10[slice]] = curY;
1865 fZ10[slice][fN10[slice]] = fZ[i];
1866 fClusterIndex10[slice][fN10[slice]]=i;
1871 for (Int_t slice=0; slice<21;slice++){
1872 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1873 fClusters20[slice][fN20[slice]] = fClusters[i];
1874 fY20[slice][fN20[slice]] = curY;
1875 fZ20[slice][fN20[slice]] = fZ[i];
1876 fClusterIndex20[slice][fN20[slice]]=i;
1883 // consistency check
1885 for (Int_t i=0;i<fN-1;i++){
1891 for (Int_t slice=0;slice<21;slice++)
1892 for (Int_t i=0;i<fN20[slice]-1;i++){
1893 if (fZ20[slice][i]>fZ20[slice][i+1]){
1900 //------------------------------------------------------------------------
1901 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1902 //--------------------------------------------------------------------
1903 // This function returns the index of the nearest cluster
1904 //--------------------------------------------------------------------
1907 if (fCurrentSlice<0) {
1916 if (ncl==0) return 0;
1917 Int_t b=0, e=ncl-1, m=(b+e)/2;
1918 for (; b<e; m=(b+e)/2) {
1919 // if (z > fClusters[m]->GetZ()) b=m+1;
1920 if (z > zcl[m]) b=m+1;
1925 //------------------------------------------------------------------------
1926 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 {
1927 //--------------------------------------------------------------------
1928 // This function computes the rectangular road for this track
1929 //--------------------------------------------------------------------
1932 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1933 // take into account the misalignment: propagate track to misaligned detector plane
1934 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1936 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1937 TMath::Sqrt(track->GetSigmaZ2() +
1938 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1939 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1940 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1941 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1942 TMath::Sqrt(track->GetSigmaY2() +
1943 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1944 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1945 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1947 // track at boundary between detectors, enlarge road
1948 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1949 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1950 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1951 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1952 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1953 Float_t tgl = TMath::Abs(track->GetTgl());
1954 if (tgl > 1.) tgl=1.;
1955 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1956 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1957 Float_t snp = TMath::Abs(track->GetSnp());
1958 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1959 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1962 // add to the road a term (up to 2-3 mm) to deal with misalignments
1963 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1964 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1966 Double_t r = fgLayers[ilayer].GetR();
1967 zmin = track->GetZ() - dz;
1968 zmax = track->GetZ() + dz;
1969 ymin = track->GetY() + r*det.GetPhi() - dy;
1970 ymax = track->GetY() + r*det.GetPhi() + dy;
1972 // bring track back to idead detector plane
1973 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1977 //------------------------------------------------------------------------
1978 void AliITStrackerMI::AliITSlayer::
1979 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1980 //--------------------------------------------------------------------
1981 // This function sets the "window"
1982 //--------------------------------------------------------------------
1984 Double_t circle=2*TMath::Pi()*fR;
1990 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1991 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1992 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1993 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1994 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1996 Float_t ymiddle = (fYmin+fYmax)*0.5;
1997 if (ymiddle<fYB[0]) {
1998 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1999 } else if (ymiddle>fYB[1]) {
2000 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
2006 fClustersCs = fClusters;
2007 fClusterIndexCs = fClusterIndex;
2013 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2014 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2015 if (slice<0) slice=0;
2016 if (slice>20) slice=20;
2017 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2019 fCurrentSlice=slice;
2020 fClustersCs = fClusters20[fCurrentSlice];
2021 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2022 fYcs = fY20[fCurrentSlice];
2023 fZcs = fZ20[fCurrentSlice];
2024 fNcs = fN20[fCurrentSlice];
2029 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2030 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2031 if (slice<0) slice=0;
2032 if (slice>10) slice=10;
2033 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2035 fCurrentSlice=slice;
2036 fClustersCs = fClusters10[fCurrentSlice];
2037 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2038 fYcs = fY10[fCurrentSlice];
2039 fZcs = fZ10[fCurrentSlice];
2040 fNcs = fN10[fCurrentSlice];
2045 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2046 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2047 if (slice<0) slice=0;
2048 if (slice>5) slice=5;
2049 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2051 fCurrentSlice=slice;
2052 fClustersCs = fClusters5[fCurrentSlice];
2053 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2054 fYcs = fY5[fCurrentSlice];
2055 fZcs = fZ5[fCurrentSlice];
2056 fNcs = fN5[fCurrentSlice];
2060 fI = FindClusterIndex(fZmin);
2061 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2067 //------------------------------------------------------------------------
2068 Int_t AliITStrackerMI::AliITSlayer::
2069 FindDetectorIndex(Double_t phi, Double_t z) const {
2070 //--------------------------------------------------------------------
2071 //This function finds the detector crossed by the track
2072 //--------------------------------------------------------------------
2074 if (fZOffset<0) // old geometry
2075 dphi = -(phi-fPhiOffset);
2076 else // new geometry
2077 dphi = phi-fPhiOffset;
2080 if (dphi < 0) dphi += 2*TMath::Pi();
2081 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2082 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2083 if (np>=fNladders) np-=fNladders;
2084 if (np<0) np+=fNladders;
2087 Double_t dz=fZOffset-z;
2088 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2089 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2090 if (nz>=fNdetectors || nz<0) {
2091 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2095 // ad hoc correction for 3rd ladder of SDD inner layer,
2096 // which is reversed (rotated by pi around local y)
2097 // this correction is OK only from AliITSv11Hybrid onwards
2098 if (GetR()>12. && GetR()<20.) { // SDD inner
2099 if(np==2) { // 3rd ladder
2100 Double_t posMod252[3];
2101 AliITSgeomTGeo::GetTranslation(252,posMod252);
2102 // check the Z coordinate of Mod 252: if negative
2103 // (old SDD geometry in AliITSv11Hybrid)
2104 // the swap of numeration whould be applied
2105 if(posMod252[2]<0.){
2106 nz = (fNdetectors-1) - nz;
2110 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2113 return np*fNdetectors + nz;
2115 //------------------------------------------------------------------------
2116 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2118 //--------------------------------------------------------------------
2119 // This function returns clusters within the "window"
2120 //--------------------------------------------------------------------
2122 if (fCurrentSlice<0) {
2123 Double_t rpi2 = 2.*fR*TMath::Pi();
2124 for (Int_t i=fI; i<fImax; i++) {
2127 if (fYmax<y) y -= rpi2;
2128 if (fYmin>y) y += rpi2;
2129 if (y<fYmin) continue;
2130 if (y>fYmax) continue;
2132 // skip clusters that are in "extended" road but they
2133 // 3sigma error does not touch the original road
2134 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2135 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2137 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2140 return fClusters[i];
2143 for (Int_t i=fI; i<fImax; i++) {
2144 if (fYcs[i]<fYmin) continue;
2145 if (fYcs[i]>fYmax) continue;
2146 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2147 ci=fClusterIndexCs[i];
2149 return fClustersCs[i];
2154 //------------------------------------------------------------------------
2155 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2157 //--------------------------------------------------------------------
2158 // This function returns the layer thickness at this point (units X0)
2159 //--------------------------------------------------------------------
2161 x0=AliITSRecoParam::GetX0Air();
2162 if (43<fR&&fR<45) { //SSD2
2165 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2166 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2167 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2168 for (Int_t i=0; i<12; i++) {
2169 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2170 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2174 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2175 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2179 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2180 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2183 if (37<fR&&fR<41) { //SSD1
2186 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2187 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2188 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2189 for (Int_t i=0; i<11; i++) {
2190 if (TMath::Abs(z-3.9*i)<0.15) {
2191 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2195 if (TMath::Abs(z+3.9*i)<0.15) {
2196 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2200 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2201 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2204 if (13<fR&&fR<26) { //SDD
2207 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2209 if (TMath::Abs(y-1.80)<0.55) {
2211 for (Int_t j=0; j<20; j++) {
2212 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2213 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2216 if (TMath::Abs(y+1.80)<0.55) {
2218 for (Int_t j=0; j<20; j++) {
2219 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2220 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2224 for (Int_t i=0; i<4; i++) {
2225 if (TMath::Abs(z-7.3*i)<0.60) {
2227 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2230 if (TMath::Abs(z+7.3*i)<0.60) {
2232 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2237 if (6<fR&&fR<8) { //SPD2
2238 Double_t dd=0.0063; x0=21.5;
2240 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2241 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2243 if (3<fR&&fR<5) { //SPD1
2244 Double_t dd=0.0063; x0=21.5;
2246 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2247 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2252 //------------------------------------------------------------------------
2253 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2255 fRmisal(det.fRmisal),
2257 fSinPhi(det.fSinPhi),
2258 fCosPhi(det.fCosPhi),
2264 fNChips(det.fNChips),
2265 fChipIsBad(det.fChipIsBad)
2269 //------------------------------------------------------------------------
2270 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2271 const AliITSDetTypeRec *detTypeRec)
2273 //--------------------------------------------------------------------
2274 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2275 //--------------------------------------------------------------------
2277 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2278 // while in the tracker they start from 0 for each layer
2279 for(Int_t il=0; il<ilayer; il++)
2280 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2283 if (ilayer==0 || ilayer==1) { // ---------- SPD
2285 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2287 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2290 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2294 // Get calibration from AliITSDetTypeRec
2295 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2296 calib->SetModuleIndex(idet);
2297 AliITSCalibration *calibSPDdead = 0;
2298 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2299 if (calib->IsBad() ||
2300 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2303 // printf("lay %d bad %d\n",ilayer,idet);
2306 // Get segmentation from AliITSDetTypeRec
2307 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2309 // Read info about bad chips
2310 fNChips = segm->GetMaximumChipIndex()+1;
2311 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2312 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2313 fChipIsBad = new Bool_t[fNChips];
2314 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2315 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2316 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2317 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2322 //------------------------------------------------------------------------
2323 Double_t AliITStrackerMI::GetEffectiveThickness()
2325 //--------------------------------------------------------------------
2326 // Returns the thickness between the current layer and the vertex (units X0)
2327 //--------------------------------------------------------------------
2330 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2331 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2332 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2336 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2337 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2341 Double_t xn=fgLayers[fI].GetR();
2342 for (Int_t i=0; i<fI; i++) {
2343 Double_t xi=fgLayers[i].GetR();
2344 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2350 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2351 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2354 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2355 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2360 //------------------------------------------------------------------------
2361 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2363 //--------------------------------------------------------------------
2364 // Returns the array of layers between the current layer and the vertex
2365 //--------------------------------------------------------------------
2368 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2369 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2370 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2375 for (int il=fI;il--;) {
2378 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2379 xMS[nl++] = AliITSRecoParam::GetrInsideShield(1);
2382 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2383 xMS[nl++] = AliITSRecoParam::GetrInsideShield(0);
2386 x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2387 xMS[nl++] = fgLayers[il].GetR();
2392 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2393 xMS[nl++] = AliITSRecoParam::GetrPipe();
2399 //------------------------------------------------------------------------
2400 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2401 //-------------------------------------------------------------------
2402 // This function returns number of clusters within the "window"
2403 //--------------------------------------------------------------------
2405 for (Int_t i=fI; i<fN; i++) {
2406 const AliITSRecPoint *c=fClusters[i];
2407 if (c->GetZ() > fZmax) break;
2408 if (c->IsUsed()) continue;
2409 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2410 Double_t y=fR*det.GetPhi() + c->GetY();
2412 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2413 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2415 if (y<fYmin) continue;
2416 if (y>fYmax) continue;
2421 //------------------------------------------------------------------------
2422 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2423 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2425 //--------------------------------------------------------------------
2426 // This function refits the track "track" at the position "x" using
2427 // the clusters from "clusters"
2428 // If "extra"==kTRUE,
2429 // the clusters from overlapped modules get attached to "track"
2430 // If "planeff"==kTRUE,
2431 // special approach for plane efficiency evaluation is applyed
2432 //--------------------------------------------------------------------
2434 Int_t index[AliITSgeomTGeo::kNLayers];
2436 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2437 Int_t nc=clusters->GetNumberOfClusters();
2438 for (k=0; k<nc; k++) {
2439 Int_t idx=clusters->GetClusterIndex(k);
2440 Int_t ilayer=(idx&0xf0000000)>>28;
2444 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2446 //------------------------------------------------------------------------
2447 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2448 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2450 //--------------------------------------------------------------------
2451 // This function refits the track "track" at the position "x" using
2452 // the clusters from array
2453 // If "extra"==kTRUE,
2454 // the clusters from overlapped modules get attached to "track"
2455 // If "planeff"==kTRUE,
2456 // special approach for plane efficiency evaluation is applyed
2457 //--------------------------------------------------------------------
2458 Int_t index[AliITSgeomTGeo::kNLayers];
2460 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2462 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2463 index[k]=clusters[k];
2466 // special for cosmics and TPC prolonged tracks:
2467 // propagate to the innermost of:
2468 // - innermost layer crossed by the track
2469 // - innermost layer where a cluster was associated to the track
2470 static AliITSRecoParam *repa = NULL;
2472 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2474 repa = AliITSRecoParam::GetHighFluxParam();
2475 AliWarning("Using default AliITSRecoParam class");
2478 Int_t evsp=repa->GetEventSpecie();
2480 if(track->GetESDtrack()) trStatus=track->GetStatus();
2481 Int_t innermostlayer=0;
2482 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2484 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2485 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2486 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2487 index[innermostlayer] >= 0 ) break;
2490 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2493 Int_t modstatus=1; // found
2495 Int_t from, to, step;
2496 if (xx > track->GetX()) {
2497 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2500 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2503 TString dir = (step>0 ? "outward" : "inward");
2505 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2506 AliITSlayer &layer=fgLayers[ilayer];
2507 Double_t r=layer.GetR();
2509 if (step<0 && xx>r) break;
2511 // material between SSD and SDD, SDD and SPD
2512 Double_t hI=ilayer-0.5*step;
2513 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2514 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2515 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2516 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2519 Double_t oldGlobXYZ[3];
2520 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2522 // continue if we are already beyond this layer
2523 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2524 if(step>0 && oldGlobR > r) continue; // going outward
2525 if(step<0 && oldGlobR < r) continue; // going inward
2528 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2530 Int_t idet=layer.FindDetectorIndex(phi,z);
2532 // check if we allow a prolongation without point for large-eta tracks
2533 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2535 modstatus = 4; // out in z
2536 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2537 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2540 // apply correction for material of the current layer
2541 // add time if going outward
2542 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2546 if (idet<0) return kFALSE;
2549 const AliITSdetector &det=layer.GetDetector(idet);
2550 // only for ITS-SA tracks refit
2551 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2553 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2555 track->SetDetectorIndex(idet);
2556 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2558 Double_t dz,zmin,zmax,dy,ymin,ymax;
2560 const AliITSRecPoint *clAcc=0;
2561 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2563 Int_t idx=index[ilayer];
2564 if (idx>=0) { // cluster in this layer
2565 modstatus = 6; // no refit
2566 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2568 if (idet != cl->GetDetectorIndex()) {
2569 idet=cl->GetDetectorIndex();
2570 const AliITSdetector &detc=layer.GetDetector(idet);
2571 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2572 track->SetDetectorIndex(idet);
2573 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2575 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2576 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2580 modstatus = 1; // found
2585 } else { // no cluster in this layer
2587 modstatus = 3; // skipped
2588 // Plane Eff determination:
2589 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2590 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2591 UseTrackForPlaneEff(track,ilayer);
2594 modstatus = 5; // no cls in road
2596 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2597 dz = 0.5*(zmax-zmin);
2598 dy = 0.5*(ymax-ymin);
2599 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2600 if (dead==1) modstatus = 7; // holes in z in SPD
2601 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2606 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2607 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2609 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2612 if (extra && clAcc) { // search for extra clusters in overlapped modules
2613 AliITStrackV2 tmp(*track);
2614 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2615 layer.SelectClusters(zmin,zmax,ymin,ymax);
2617 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2619 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2620 Double_t tolerance=0.1;
2621 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2622 // only clusters in another module! (overlaps)
2623 idetExtra = clExtra->GetDetectorIndex();
2624 if (idet == idetExtra) continue;
2626 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2628 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2629 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2630 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2631 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2633 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2634 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2637 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2638 track->SetExtraModule(ilayer,idetExtra);
2640 } // end search for extra clusters in overlapped modules
2642 // Correct for material of the current layer
2644 // add time if going outward
2645 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2646 track->SetCheckInvariant(kTRUE);
2647 } // end loop on layers
2649 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2653 //------------------------------------------------------------------------
2654 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2657 // calculate normalized chi2
2658 // return NormalizedChi2(track,0);
2661 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2662 // track->fdEdxMismatch=0;
2663 Float_t dedxmismatch =0;
2664 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2666 for (Int_t i = 0;i<6;i++){
2667 if (track->GetClIndex(i)>=0){
2668 Float_t cerry, cerrz;
2669 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2671 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2674 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2675 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2676 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2678 cchi2+=(0.5-ratio)*10.;
2679 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2680 dedxmismatch+=(0.5-ratio)*10.;
2684 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2685 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2686 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2687 if (i<2) chi2+=2*cl->GetDeltaProbability();
2693 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2694 track->SetdEdxMismatch(dedxmismatch);
2698 for (Int_t i = 0;i<4;i++){
2699 if (track->GetClIndex(i)>=0){
2700 Float_t cerry, cerrz;
2701 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2702 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2705 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2706 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2710 for (Int_t i = 4;i<6;i++){
2711 if (track->GetClIndex(i)>=0){
2712 Float_t cerry, cerrz;
2713 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2714 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2717 Float_t cerryb, cerrzb;
2718 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2719 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2722 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2723 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2728 if (track->GetESDtrack()->GetTPCsignal()>85){
2729 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2731 chi2+=(0.5-ratio)*5.;
2734 chi2+=(ratio-2.0)*3;
2738 Double_t match = TMath::Sqrt(track->GetChi22());
2739 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2740 if (!track->GetConstrain()) {
2741 if (track->GetNumberOfClusters()>2) {
2742 match/=track->GetNumberOfClusters()-2.;
2747 if (match<0) match=0;
2749 // penalty factor for missing points (NDeadZone>0), but no penalty
2750 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2751 // or there is a dead from OCDB)
2752 Float_t deadzonefactor = 0.;
2753 if (track->GetNDeadZone()>0.) {
2754 Int_t sumDeadZoneProbability=0;
2755 for(Int_t ilay=0;ilay<6;ilay++) {
2756 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2758 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2759 if(nDeadZoneWithProbNot1>0) {
2760 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2761 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2762 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2764 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2765 deadzonefactor = 3.*(1.1-deadZoneProbability);
2769 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2770 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2771 1./(1.+track->GetNSkipped()));
2772 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2773 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2776 //------------------------------------------------------------------------
2777 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2780 // return matching chi2 between two tracks
2781 Double_t largeChi2=1000.;
2783 AliITStrackMI track3(*track2);
2784 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2786 vec(0,0)=track1->GetY() - track3.GetY();
2787 vec(1,0)=track1->GetZ() - track3.GetZ();
2788 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2789 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2790 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2793 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2794 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2795 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2796 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2797 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2799 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2800 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2801 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2802 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2804 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2805 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2806 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2808 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2809 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2811 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2814 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2815 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2818 //------------------------------------------------------------------------
2819 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2822 // return probability that given point (characterized by z position and error)
2823 // is in SPD dead zone
2824 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2826 Double_t probability = 0.;
2827 Double_t nearestz = 0.,distz=0.;
2828 Int_t nearestzone = -1;
2829 Double_t mindistz = 1000.;
2831 // find closest dead zone
2832 for (Int_t i=0; i<3; i++) {
2833 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2834 if (distz<mindistz) {
2836 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2841 // too far from dead zone
2842 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2845 Double_t zmin, zmax;
2846 if (nearestzone==0) { // dead zone at z = -7
2847 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2848 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2849 } else if (nearestzone==1) { // dead zone at z = 0
2850 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2851 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2852 } else if (nearestzone==2) { // dead zone at z = +7
2853 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2854 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2859 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2861 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2862 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2863 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2866 //------------------------------------------------------------------------
2867 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2870 // calculate normalized chi2
2872 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2874 for (Int_t i = 0;i<6;i++){
2875 if (TMath::Abs(track->GetDy(i))>0){
2876 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2877 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2880 else{chi2[i]=10000;}
2883 TMath::Sort(6,chi2,index,kFALSE);
2884 Float_t max = float(ncl)*fac-1.;
2885 Float_t sumchi=0, sumweight=0;
2886 for (Int_t i=0;i<max+1;i++){
2887 Float_t weight = (i<max)?1.:(max+1.-i);
2888 sumchi+=weight*chi2[index[i]];
2891 Double_t normchi2 = sumchi/sumweight;
2894 //------------------------------------------------------------------------
2895 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2898 // calculate normalized chi2
2899 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2902 for (Int_t i=0;i<6;i++){
2903 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2904 Double_t sy1 = forwardtrack->GetSigmaY(i);
2905 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2906 Double_t sy2 = backtrack->GetSigmaY(i);
2907 Double_t sz2 = backtrack->GetSigmaZ(i);
2908 if (i<2){ sy2=1000.;sz2=1000;}
2910 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2911 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2913 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2914 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2916 res+= nz0*nz0+ny0*ny0;
2919 if (npoints>1) return
2920 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2921 //2*forwardtrack->fNUsed+
2922 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2923 1./(1.+forwardtrack->GetNSkipped()));
2926 //------------------------------------------------------------------------
2927 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2928 //--------------------------------------------------------------------
2929 // Return pointer to a given cluster
2930 //--------------------------------------------------------------------
2931 Int_t l=(index & 0xf0000000) >> 28;
2932 Int_t c=(index & 0x0fffffff) >> 00;
2933 return fgLayers[l].GetWeight(c);
2935 //------------------------------------------------------------------------
2936 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2938 //---------------------------------------------
2939 // register track to the list
2941 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2944 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2945 Int_t index = track->GetClusterIndex(icluster);
2946 Int_t l=(index & 0xf0000000) >> 28;
2947 Int_t c=(index & 0x0fffffff) >> 00;
2948 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2949 for (Int_t itrack=0;itrack<4;itrack++){
2950 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2951 fgLayers[l].SetClusterTracks(itrack,c,id);
2957 //------------------------------------------------------------------------
2958 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2960 //---------------------------------------------
2961 // unregister track from the list
2962 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2963 Int_t index = track->GetClusterIndex(icluster);
2964 Int_t l=(index & 0xf0000000) >> 28;
2965 Int_t c=(index & 0x0fffffff) >> 00;
2966 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2967 for (Int_t itrack=0;itrack<4;itrack++){
2968 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2969 fgLayers[l].SetClusterTracks(itrack,c,-1);
2974 //------------------------------------------------------------------------
2975 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2977 //-------------------------------------------------------------
2978 //get number of shared clusters
2979 //-------------------------------------------------------------
2981 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2982 // mean number of clusters
2983 Float_t *ny = GetNy(id), *nz = GetNz(id);
2986 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2987 Int_t index = track->GetClusterIndex(icluster);
2988 Int_t l=(index & 0xf0000000) >> 28;
2989 Int_t c=(index & 0x0fffffff) >> 00;
2990 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2991 // if (ny[l]<1.e-13){
2992 // printf("problem\n");
2994 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2998 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2999 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3000 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3002 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3005 deltan = (cl->GetNz()-nz[l]);
3007 if (deltan>2.0) continue; // extended - highly probable shared cluster
3008 weight = 2./TMath::Max(3.+deltan,2.);
3010 for (Int_t itrack=0;itrack<4;itrack++){
3011 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3013 clist[l] = (AliITSRecPoint*)GetCluster(index);
3014 track->SetSharedWeight(l,weight);
3020 track->SetNUsed(shared);
3023 //------------------------------------------------------------------------
3024 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3027 // find first shared track
3029 // mean number of clusters
3030 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
3032 for (Int_t i=0;i<6;i++) overlist[i]=-1;
3033 Int_t sharedtrack=100000;
3034 Int_t tracks[24],trackindex=0;
3035 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3037 for (Int_t icluster=0;icluster<6;icluster++){
3038 if (clusterlist[icluster]<0) continue;
3039 Int_t index = clusterlist[icluster];
3040 Int_t l=(index & 0xf0000000) >> 28;
3041 Int_t c=(index & 0x0fffffff) >> 00;
3042 // if (ny[l]<1.e-13){
3043 // printf("problem\n");
3045 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3046 //if (l>3) continue;
3047 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3050 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3051 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3052 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3054 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3057 deltan = (cl->GetNz()-nz[l]);
3059 if (deltan>2.0) continue; // extended - highly probable shared cluster
3061 for (Int_t itrack=3;itrack>=0;itrack--){
3062 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3063 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3064 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
3069 if (trackindex==0) return -1;
3071 sharedtrack = tracks[0];
3073 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3076 Int_t tracks2[24], cluster[24];
3077 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
3080 for (Int_t i=0;i<trackindex;i++){
3081 if (tracks[i]<0) continue;
3082 tracks2[index] = tracks[i];
3084 for (Int_t j=i+1;j<trackindex;j++){
3085 if (tracks[j]<0) continue;
3086 if (tracks[j]==tracks[i]){
3094 for (Int_t i=0;i<index;i++){
3095 if (cluster[index]>max) {
3096 sharedtrack=tracks2[index];
3103 if (sharedtrack>=100000) return -1;
3105 // make list of overlaps
3107 for (Int_t icluster=0;icluster<6;icluster++){
3108 if (clusterlist[icluster]<0) continue;
3109 Int_t index = clusterlist[icluster];
3110 Int_t l=(index & 0xf0000000) >> 28;
3111 Int_t c=(index & 0x0fffffff) >> 00;
3112 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3113 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3115 if (cl->GetNy()>2) continue;
3116 if (cl->GetNz()>2) continue;
3119 if (cl->GetNy()>3) continue;
3120 if (cl->GetNz()>3) continue;
3123 for (Int_t itrack=3;itrack>=0;itrack--){
3124 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3125 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3133 //------------------------------------------------------------------------
3134 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3136 // try to find track hypothesys without conflicts
3137 // with minimal chi2;
3138 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3139 Int_t entries1 = arr1->GetEntriesFast();
3140 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3141 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3142 Int_t entries2 = arr2->GetEntriesFast();
3143 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3145 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3146 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3147 if (track10->Pt()>0.5+track20->Pt()) return track10;
3149 for (Int_t itrack=0;itrack<entries1;itrack++){
3150 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3151 UnRegisterClusterTracks(track,trackID1);
3154 for (Int_t itrack=0;itrack<entries2;itrack++){
3155 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3156 UnRegisterClusterTracks(track,trackID2);
3160 Float_t maxconflicts=6;
3161 Double_t maxchi2 =1000.;
3163 // get weight of hypothesys - using best hypothesy
3166 Int_t list1[6],list2[6];
3167 AliITSRecPoint *clist1[6], *clist2[6] ;
3168 RegisterClusterTracks(track10,trackID1);
3169 RegisterClusterTracks(track20,trackID2);
3170 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3171 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3172 UnRegisterClusterTracks(track10,trackID1);
3173 UnRegisterClusterTracks(track20,trackID2);
3176 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3177 Float_t nerry[6],nerrz[6];
3178 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3179 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3180 for (Int_t i=0;i<6;i++){
3181 if ( (erry1[i]>0) && (erry2[i]>0)) {
3182 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3183 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3185 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3186 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3188 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3189 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3190 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3193 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3194 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3195 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3203 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3204 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3205 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3206 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3208 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3209 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3210 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3212 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3213 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3214 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3217 Double_t sumw = w1+w2;
3221 w1 = (d2+0.5)/(d1+d2+1.);
3222 w2 = (d1+0.5)/(d1+d2+1.);
3224 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3225 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3227 // get pair of "best" hypothesys
3229 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3230 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3232 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3233 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3234 //if (track1->fFakeRatio>0) continue;
3235 RegisterClusterTracks(track1,trackID1);
3236 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3237 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3239 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3240 //if (track2->fFakeRatio>0) continue;
3242 RegisterClusterTracks(track2,trackID2);
3243 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3244 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3245 UnRegisterClusterTracks(track2,trackID2);
3247 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3248 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3249 if (nskipped>0.5) continue;
3251 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3252 if (conflict1+1<cconflict1) continue;
3253 if (conflict2+1<cconflict2) continue;
3257 for (Int_t i=0;i<6;i++){
3259 Float_t c1 =0.; // conflict coeficients
3261 if (clist1[i]&&clist2[i]){
3264 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3267 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3269 c1 = 2./TMath::Max(3.+deltan,2.);
3270 c2 = 2./TMath::Max(3.+deltan,2.);
3276 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3279 deltan = (clist1[i]->GetNz()-nz1[i]);
3281 c1 = 2./TMath::Max(3.+deltan,2.);
3288 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3291 deltan = (clist2[i]->GetNz()-nz2[i]);
3293 c2 = 2./TMath::Max(3.+deltan,2.);
3299 if (TMath::Abs(track1->GetDy(i))>0.) {
3300 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3301 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3302 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3303 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3305 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3308 if (TMath::Abs(track2->GetDy(i))>0.) {
3309 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3310 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3311 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3312 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3315 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3317 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3318 if (chi21>0) sum+=w1;
3319 if (chi22>0) sum+=w2;
3322 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3323 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3324 Double_t normchi2 = 2*conflict+sumchi2/sum;
3325 if ( normchi2 <maxchi2 ){
3328 maxconflicts = conflict;
3332 UnRegisterClusterTracks(track1,trackID1);
3335 // if (maxconflicts<4 && maxchi2<th0){
3336 if (maxchi2<th0*2.){
3337 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3338 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3339 track1->SetChi2MIP(5,maxconflicts);
3340 track1->SetChi2MIP(6,maxchi2);
3341 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3342 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3343 track1->SetChi2MIP(8,index1);
3344 fBestTrackIndex[trackID1] =index1;
3345 UpdateESDtrack(track1, AliESDtrack::kITSin);
3346 original->SetWinner(track1);
3348 else if (track10->GetChi2MIP(0)<th1){
3349 track10->SetChi2MIP(5,maxconflicts);
3350 track10->SetChi2MIP(6,maxchi2);
3351 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3352 UpdateESDtrack(track10,AliESDtrack::kITSin);
3353 original->SetWinner(track10);
3356 for (Int_t itrack=0;itrack<entries1;itrack++){
3357 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3358 UnRegisterClusterTracks(track,trackID1);
3361 for (Int_t itrack=0;itrack<entries2;itrack++){
3362 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3363 UnRegisterClusterTracks(track,trackID2);
3366 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3367 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3368 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3369 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3370 RegisterClusterTracks(track10,trackID1);
3372 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3373 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3374 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3375 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3376 RegisterClusterTracks(track20,trackID2);
3381 //------------------------------------------------------------------------
3382 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3383 //--------------------------------------------------------------------
3384 // This function marks clusters assigned to the track
3385 //--------------------------------------------------------------------
3386 AliTracker::UseClusters(t,from);
3388 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3389 //if (c->GetQ()>2) c->Use();
3390 if (c->GetSigmaZ2()>0.1) c->Use();
3391 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3392 //if (c->GetQ()>2) c->Use();
3393 if (c->GetSigmaZ2()>0.1) c->Use();
3396 //------------------------------------------------------------------------
3397 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3399 //------------------------------------------------------------------
3400 // add track to the list of hypothesys
3401 //------------------------------------------------------------------
3403 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3404 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3406 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3408 array = new TObjArray(10);
3409 fTrackHypothesys.AddAt(array,esdindex);
3411 array->AddLast(track);
3413 //------------------------------------------------------------------------
3414 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3416 //-------------------------------------------------------------------
3417 // compress array of track hypothesys
3418 // keep only maxsize best hypothesys
3419 //-------------------------------------------------------------------
3420 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3421 if (! (fTrackHypothesys.At(esdindex)) ) return;
3422 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3423 Int_t entries = array->GetEntriesFast();
3425 //- find preliminary besttrack as a reference
3426 Float_t minchi2=10000;
3428 AliITStrackMI * besttrack=0;
3430 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3431 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3432 if (!track) continue;
3433 Float_t chi2 = NormalizedChi2(track,0);
3435 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3436 track->SetLabel(tpcLabel);
3438 track->SetFakeRatio(1.);
3439 CookLabel(track,0.); //For comparison only
3441 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3442 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3443 if (track->GetNumberOfClusters()<maxn) continue;
3444 maxn = track->GetNumberOfClusters();
3445 // if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3452 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3453 delete array->RemoveAt(itrack);
3457 if (!besttrack) return;
3460 //take errors of best track as a reference
3461 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3462 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3463 for (Int_t j=0;j<6;j++) {
3464 if (besttrack->GetClIndex(j)>=0){
3465 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3466 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3467 ny[j] = besttrack->GetNy(j);
3468 nz[j] = besttrack->GetNz(j);
3472 // calculate normalized chi2
3474 Float_t * chi2 = new Float_t[entries];
3475 Int_t * index = new Int_t[entries];
3476 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3477 for (Int_t itrack=0;itrack<entries;itrack++){
3478 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3480 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3481 double chi2t = GetNormalizedChi2(track, mode);
3482 track->SetChi2MIP(0,chi2t);
3483 if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3484 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3485 chi2[itrack] = chi2t;
3488 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3489 delete array->RemoveAt(itrack);
3495 TMath::Sort(entries,chi2,index,kFALSE);
3496 besttrack = (AliITStrackMI*)array->At(index[0]);
3497 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3498 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3499 for (Int_t j=0;j<6;j++){
3500 if (besttrack->GetClIndex(j)>=0){
3501 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3502 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3503 ny[j] = besttrack->GetNy(j);
3504 nz[j] = besttrack->GetNz(j);
3509 // calculate one more time with updated normalized errors
3510 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3511 for (Int_t itrack=0;itrack<entries;itrack++){
3512 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3514 double chi2t = GetNormalizedChi2(track, mode);
3515 track->SetChi2MIP(0,chi2t);
3516 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3517 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3518 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3519 chi2[itrack] = chi2t; //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3522 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3523 delete array->RemoveAt(itrack);
3528 entries = array->GetEntriesFast();
3532 TObjArray * newarray = new TObjArray();
3533 TMath::Sort(entries,chi2,index,kFALSE);
3534 besttrack = (AliITStrackMI*)array->At(index[0]);
3536 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3538 for (Int_t j=0;j<6;j++){
3539 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3540 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3541 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3542 ny[j] = besttrack->GetNy(j);
3543 nz[j] = besttrack->GetNz(j);
3546 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3547 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3548 Float_t minn = besttrack->GetNumberOfClusters()-3;
3550 for (Int_t i=0;i<entries;i++){
3551 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3552 if (!track) continue;
3553 if (accepted>maxcut) break;
3554 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3555 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3556 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3557 delete array->RemoveAt(index[i]);
3561 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3562 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3563 if (!shortbest) accepted++;
3565 newarray->AddLast(array->RemoveAt(index[i]));
3566 for (Int_t j=0;j<6;j++){
3568 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3569 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3570 ny[j] = track->GetNy(j);
3571 nz[j] = track->GetNz(j);
3576 delete array->RemoveAt(index[i]);
3580 delete fTrackHypothesys.RemoveAt(esdindex);
3581 fTrackHypothesys.AddAt(newarray,esdindex);
3585 delete fTrackHypothesys.RemoveAt(esdindex);
3591 //------------------------------------------------------------------------
3592 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3594 //-------------------------------------------------------------
3595 // try to find best hypothesy
3596 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3597 // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3598 //-------------------------------------------------------------
3599 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3600 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3601 if (!array) return 0;
3602 Int_t entries = array->GetEntriesFast();
3603 if (!entries) return 0;
3604 Float_t minchi2 = 100000;
3605 AliITStrackMI * besttrack=0;
3607 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3608 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3609 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3610 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3612 for (Int_t i=0;i<entries;i++){
3613 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3614 if (!track) continue;
3615 Float_t sigmarfi,sigmaz;
3616 GetDCASigma(track,sigmarfi,sigmaz);
3617 track->SetDnorm(0,sigmarfi);
3618 track->SetDnorm(1,sigmaz);
3620 track->SetChi2MIP(1,1000000);
3621 track->SetChi2MIP(2,1000000);
3622 track->SetChi2MIP(3,1000000);
3625 backtrack = new(backtrack) AliITStrackMI(*track);
3626 if (track->GetConstrain()) {
3627 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3628 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3629 if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3630 else {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3632 backtrack->ResetCovariance(10.);
3634 backtrack->ResetCovariance(10.);
3636 backtrack->ResetClusters();
3638 Double_t x = original->GetX();
3639 if (!RefitAt(x,backtrack,track)) continue;
3641 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3642 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3643 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3644 track->SetChi22(GetMatchingChi2(backtrack,original));
3646 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3647 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3648 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3651 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3653 //forward track - without constraint
3654 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3655 forwardtrack->ResetClusters();
3657 if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue; // w/o fwd track MIP03 is meaningless
3658 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3659 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3660 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3662 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3663 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3664 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3665 forwardtrack->SetD(0,track->GetD(0));
3666 forwardtrack->SetD(1,track->GetD(1));
3669 AliITSRecPoint* clist[6];
3670 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3671 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3674 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3675 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3676 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3677 track->SetChi2MIP(3,1000);
3680 Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed(); //RS
3681 if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3682 else chi2 += track->GetNUsed();
3684 for (Int_t ichi=0;ichi<5;ichi++){
3685 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3687 if (chi2 < minchi2){
3688 //besttrack = new AliITStrackMI(*forwardtrack);
3690 besttrack->SetLabel(track->GetLabel());
3691 besttrack->SetFakeRatio(track->GetFakeRatio());
3693 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3694 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3695 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3699 delete forwardtrack;
3701 if (!besttrack) return 0;
3704 for (Int_t i=0;i<entries;i++){
3705 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3707 if (!track) continue;
3709 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3710 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3711 // RS: don't apply this cut when fSelectBestMIP03 is on
3712 || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3714 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3715 delete array->RemoveAt(i);
3725 SortTrackHypothesys(esdindex,checkmax,1);
3727 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3728 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3729 besttrack = (AliITStrackMI*)array->At(0);
3730 if (!besttrack) return 0;
3731 besttrack->SetChi2MIP(8,0);
3732 fBestTrackIndex[esdindex]=0;
3733 entries = array->GetEntriesFast();
3734 AliITStrackMI *longtrack =0;
3736 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3737 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3738 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3739 if (!track->GetConstrain()) continue;
3740 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3741 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3742 if (track->GetChi2MIP(0)>4.) continue;
3743 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3746 //if (longtrack) besttrack=longtrack;
3748 // RS do shared cluster analysis here only if the new sharing analysis is not requested
3749 //RRR if (fFlagFakes) return besttrack;
3752 AliITSRecPoint * clist[6];
3753 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3754 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3755 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3756 RegisterClusterTracks(besttrack,esdindex);
3763 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3764 if (sharedtrack>=0){
3766 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3768 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3774 if (shared>2.5) return 0;
3775 if (shared>1.0) return besttrack;
3777 // Don't sign clusters if not gold track
3779 if (!besttrack->IsGoldPrimary()) return besttrack;
3780 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3782 if (fConstraint[fPass]){
3786 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3787 for (Int_t i=0;i<6;i++){
3788 Int_t index = besttrack->GetClIndex(i);
3789 if (index<0) continue;
3790 Int_t ilayer = (index & 0xf0000000) >> 28;
3791 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3792 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3794 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3795 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3796 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3797 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3798 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3799 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3801 Bool_t cansign = kTRUE;
3802 for (Int_t itrack=0;itrack<entries; itrack++){
3803 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3804 if (!track) continue;
3805 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3806 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3812 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3813 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3814 if (!c->IsUsed()) c->Use();
3820 //------------------------------------------------------------------------
3821 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3824 // get "best" hypothesys
3827 Int_t nentries = itsTracks.GetEntriesFast();
3828 for (Int_t i=0;i<nentries;i++){
3829 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3830 if (!track) continue;
3831 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3832 if (!array) continue;
3833 if (array->GetEntriesFast()<=0) continue;
3835 AliITStrackMI* longtrack=0;
3837 Float_t maxchi2=1000;
3838 for (Int_t j=0;j<array->GetEntriesFast();j++){
3839 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3840 if (!trackHyp) continue;
3841 if (trackHyp->GetGoldV0()) {
3842 longtrack = trackHyp; //gold V0 track taken
3845 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3846 Float_t chi2 = trackHyp->GetChi2MIP(0);
3847 if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3848 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3850 if (fAfterV0){ // ??? RS
3851 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3853 if (chi2 > maxchi2) continue;
3854 minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3855 if (fSelectBestMIP03) minn++; // allow next to longest to win
3862 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3863 if (!longtrack) {longtrack = besttrack;}
3864 else besttrack= longtrack;
3868 AliITSRecPoint * clist[6];
3869 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3871 track->SetNUsed(shared);
3872 track->SetNSkipped(besttrack->GetNSkipped());
3873 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3875 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3879 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3880 //if (sharedtrack==-1) sharedtrack=0;
3881 if (sharedtrack>=0) {
3882 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3885 if (besttrack&&fAfterV0) {
3886 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3887 track->SetWinner(besttrack);
3890 if (fConstraint[fPass]) {
3891 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3892 track->SetWinner(besttrack);
3894 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3895 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3896 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3903 //------------------------------------------------------------------------
3904 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3907 // RS: flag those tracks which are suxpected to have fake clusters
3909 const double kThreshPt = 0.5;
3910 AliRefArray *refArr[6];
3912 for (int i=0;i<6;i++) {
3913 int ncl = fgLayers[i].GetNumberOfClusters();
3914 refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3916 Int_t nentries = itsTracks.GetEntriesFast();
3918 // fill cluster->track associations
3919 for (Int_t itr=0;itr<nentries;itr++){
3920 AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3921 if (!track) continue;
3922 AliITStrackMI* trackITS = track->GetWinner();
3923 if (!trackITS) continue;
3924 for (int il=trackITS->GetNumberOfClusters();il--;) {
3925 int idx = trackITS->GetClusterIndex(il);
3926 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3927 // if (c>fgLayers[l].GetNumberOfClusters()) continue;
3928 refArr[l]->AddReference(c, itr);
3932 const UInt_t kMaxRef = 100;
3933 UInt_t crefs[kMaxRef];
3935 // process tracks with shared clusters
3936 for (int itr=0;itr<nentries;itr++){
3937 AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3938 AliITStrackMI* trackH0 = track0->GetWinner();
3939 if (!trackH0) continue;
3940 AliESDtrack* esd0 = track0->GetESDtrack();
3942 for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3943 int idx = trackH0->GetClusterIndex(il);
3944 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3945 ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3946 if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3947 esd0->SetITSSharedFlag(l);
3948 for (int ir=ncrefs;ir--;) {
3949 if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3950 AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3951 AliITStrackMI* trackH1 = track1->GetWinner();
3952 AliESDtrack* esd1 = track1->GetESDtrack();
3953 esd1->SetITSSharedFlag(l);
3955 double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3956 if (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3957 else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3959 // select the one with smallest chi2's product
3960 res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3961 res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3963 if (res<0) esd1->SetITSFakeFlag(); // esd0 is winner
3964 else esd0->SetITSFakeFlag(); // esd1 is winner
3971 for (int i=6;i--;) delete refArr[i];
3974 //------------------------------------------------------------------------
3975 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3976 //--------------------------------------------------------------------
3977 //This function "cooks" a track label. If label<0, this track is fake.
3978 //--------------------------------------------------------------------
3981 if (track->GetESDtrack()){
3982 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3983 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3984 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3986 track->SetChi2MIP(9,0);
3988 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3989 Int_t cindex = track->GetClusterIndex(i);
3990 Int_t l=(cindex & 0xf0000000) >> 28;
3991 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3993 for (Int_t ind=0;ind<3;ind++){
3994 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3995 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3997 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4000 Int_t nclusters = track->GetNumberOfClusters();
4001 if (nclusters > 0) //PH Some tracks don't have any cluster
4002 track->SetFakeRatio(double(nwrong)/double(nclusters));
4003 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4004 track->SetLabel(-tpcLabel);
4006 track->SetLabel(tpcLabel);
4008 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4011 //------------------------------------------------------------------------
4012 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4014 // Fill the dE/dx in this track
4016 track->SetChi2MIP(9,0);
4017 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4018 Int_t cindex = track->GetClusterIndex(i);
4019 Int_t l=(cindex & 0xf0000000) >> 28;
4020 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4021 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4023 for (Int_t ind=0;ind<3;ind++){
4024 if (cl->GetLabel(ind)==lab) isWrong=0;
4026 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4030 track->CookdEdx(low,up);
4032 //------------------------------------------------------------------------
4033 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4035 // Create some arrays
4037 if (fCoefficients) delete []fCoefficients;
4038 fCoefficients = new Float_t[ntracks*48];
4039 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4041 //------------------------------------------------------------------------
4042 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4045 // Compute predicted chi2
4047 // Take into account the mis-alignment (bring track to cluster plane)
4048 Double_t xTrOrig=track->GetX();
4049 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4050 Float_t erry,errz,covyz;
4051 Float_t theta = track->GetTgl();
4052 Float_t phi = track->GetSnp();
4053 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4054 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4055 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()));
4056 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()));
4057 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4058 // Bring the track back to detector plane in ideal geometry
4059 // [mis-alignment will be accounted for in UpdateMI()]
4060 if (!track->Propagate(xTrOrig)) return 1000.;
4062 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4063 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4065 chi2+=0.5*TMath::Min(delta/2,2.);
4066 chi2+=2.*cluster->GetDeltaProbability();
4069 track->SetNy(layer,ny);
4070 track->SetNz(layer,nz);
4071 track->SetSigmaY(layer,erry);
4072 track->SetSigmaZ(layer, errz);
4073 track->SetSigmaYZ(layer,covyz);
4074 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4075 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4079 //------------------------------------------------------------------------
4080 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4085 Int_t layer = (index & 0xf0000000) >> 28;
4086 track->SetClIndex(layer, index);
4087 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4088 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4089 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4090 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4094 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
4097 // Take into account the mis-alignment (bring track to cluster plane)
4098 Double_t xTrOrig=track->GetX();
4099 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4100 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4101 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4102 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
4104 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4107 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4108 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4109 c.SetSigmaYZ(track->GetSigmaYZ(layer));
4112 Int_t updated = track->UpdateMI(&c,chi2,index);
4113 // Bring the track back to detector plane in ideal geometry
4114 if (!track->Propagate(xTrOrig)) return 0;
4116 if(!updated) AliDebug(2,"update failed");
4120 //------------------------------------------------------------------------
4121 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4124 //DCA sigmas parameterization
4125 //to be paramterized using external parameters in future
4128 Double_t curv=track->GetC();
4129 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4130 sigmaz = 0.0110+4.37*TMath::Abs(curv);
4132 //------------------------------------------------------------------------
4133 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4136 // Clusters from delta electrons?
4138 Int_t entries = clusterArray->GetEntriesFast();
4139 if (entries<4) return;
4140 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4141 Int_t layer = cluster->GetLayer();
4142 if (layer>1) return;
4144 Int_t ncandidates=0;
4145 Float_t r = (layer>0)? 7:4;
4147 for (Int_t i=0;i<entries;i++){
4148 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4149 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4150 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4151 index[ncandidates] = i; //candidate to belong to delta electron track
4153 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4154 cl0->SetDeltaProbability(1);
4160 for (Int_t i=0;i<ncandidates;i++){
4161 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4162 if (cl0->GetDeltaProbability()>0.8) continue;
4165 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4166 sumy=sumz=sumy2=sumyz=sumw=0.0;
4167 for (Int_t j=0;j<ncandidates;j++){
4169 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4171 Float_t dz = cl0->GetZ()-cl1->GetZ();
4172 Float_t dy = cl0->GetY()-cl1->GetY();
4173 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4174 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4175 y[ncl] = cl1->GetY();
4176 z[ncl] = cl1->GetZ();
4177 sumy+= y[ncl]*weight;
4178 sumz+= z[ncl]*weight;
4179 sumy2+=y[ncl]*y[ncl]*weight;
4180 sumyz+=y[ncl]*z[ncl]*weight;
4185 if (ncl<4) continue;
4186 Float_t det = sumw*sumy2 - sumy*sumy;
4188 if (TMath::Abs(det)>0.01){
4189 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
4190 Float_t k = (sumyz*sumw - sumy*sumz)/det;
4191 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4194 Float_t z0 = sumyz/sumy;
4195 delta = TMath::Abs(cl0->GetZ()-z0);
4198 cl0->SetDeltaProbability(1-20.*delta);
4202 //------------------------------------------------------------------------
4203 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4208 track->UpdateESDtrack(flags);
4209 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4210 if (oldtrack) delete oldtrack;
4211 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4212 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4213 // printf("Problem\n");
4216 //------------------------------------------------------------------------
4217 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4219 // Get nearest upper layer close to the point xr.
4220 // rough approximation
4222 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4223 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4225 for (Int_t i=0;i<6;i++){
4226 if (radius<kRadiuses[i]){
4233 //------------------------------------------------------------------------
4234 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4235 //--------------------------------------------------------------------
4236 // Fill a look-up table with mean material
4237 //--------------------------------------------------------------------
4241 Double_t point1[3],point2[3];
4242 Double_t phi,cosphi,sinphi,z;
4243 // 0-5 layers, 6 pipe, 7-8 shields
4244 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4245 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4247 Int_t ifirst=0,ilast=0;
4248 if(material.Contains("Pipe")) {
4250 } else if(material.Contains("Shields")) {
4252 } else if(material.Contains("Layers")) {
4255 Error("BuildMaterialLUT","Wrong layer name\n");
4258 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4259 Double_t param[5]={0.,0.,0.,0.,0.};
4260 for (Int_t i=0; i<n; i++) {
4261 phi = 2.*TMath::Pi()*gRandom->Rndm();
4262 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4263 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4264 point1[0] = rmin[imat]*cosphi;
4265 point1[1] = rmin[imat]*sinphi;
4267 point2[0] = rmax[imat]*cosphi;
4268 point2[1] = rmax[imat]*sinphi;
4270 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4271 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4273 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4275 fxOverX0Layer[imat] = param[1];
4276 fxTimesRhoLayer[imat] = param[0]*param[4];
4277 } else if(imat==6) {
4278 fxOverX0Pipe = param[1];
4279 fxTimesRhoPipe = param[0]*param[4];
4280 } else if(imat==7) {
4281 fxOverX0Shield[0] = param[1];
4282 fxTimesRhoShield[0] = param[0]*param[4];
4283 } else if(imat==8) {
4284 fxOverX0Shield[1] = param[1];
4285 fxTimesRhoShield[1] = param[0]*param[4];
4289 printf("%s\n",material.Data());
4290 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4291 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4292 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4293 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4294 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4295 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4296 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4297 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4298 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4302 //------------------------------------------------------------------------
4303 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4304 TString direction) {
4305 //-------------------------------------------------------------------
4306 // Propagate beyond beam pipe and correct for material
4307 // (material budget in different ways according to fUseTGeo value)
4308 // Add time if going outward (PropagateTo or PropagateToTGeo)
4309 //-------------------------------------------------------------------
4311 // Define budget mode:
4312 // 0: material from AliITSRecoParam (hard coded)
4313 // 1: material from TGeo in one step (on the fly)
4314 // 2: material from lut
4315 // 3: material from TGeo in one step (same for all hypotheses)
4328 if(fTrackingPhase.Contains("Clusters2Tracks"))
4329 { mode=3; } else { mode=1; }
4332 if(fTrackingPhase.Contains("Clusters2Tracks"))
4333 { mode=3; } else { mode=2; }
4339 if(fTrackingPhase.Contains("Default")) mode=0;
4341 Int_t index=fCurrentEsdTrack;
4343 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4344 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4346 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4348 Double_t xOverX0,x0,lengthTimesMeanDensity;
4352 xOverX0 = AliITSRecoParam::GetdPipe();
4353 x0 = AliITSRecoParam::GetX0Be();
4354 lengthTimesMeanDensity = xOverX0*x0;
4355 lengthTimesMeanDensity *= dir;
4356 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4359 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4362 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4363 xOverX0 = fxOverX0Pipe;
4364 lengthTimesMeanDensity = fxTimesRhoPipe;
4365 lengthTimesMeanDensity *= dir;
4366 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4369 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4370 if(fxOverX0PipeTrks[index]<0) {
4371 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4372 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4373 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4374 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4375 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4378 xOverX0 = fxOverX0PipeTrks[index];
4379 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4380 lengthTimesMeanDensity *= dir;
4381 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4387 //------------------------------------------------------------------------
4388 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4390 TString direction) {
4391 //-------------------------------------------------------------------
4392 // Propagate beyond SPD or SDD shield and correct for material
4393 // (material budget in different ways according to fUseTGeo value)
4394 // Add time if going outward (PropagateTo or PropagateToTGeo)
4395 //-------------------------------------------------------------------
4397 // Define budget mode:
4398 // 0: material from AliITSRecoParam (hard coded)
4399 // 1: material from TGeo in steps of X cm (on the fly)
4400 // X = AliITSRecoParam::GetStepSizeTGeo()
4401 // 2: material from lut
4402 // 3: material from TGeo in one step (same for all hypotheses)
4415 if(fTrackingPhase.Contains("Clusters2Tracks"))
4416 { mode=3; } else { mode=1; }
4419 if(fTrackingPhase.Contains("Clusters2Tracks"))
4420 { mode=3; } else { mode=2; }
4426 if(fTrackingPhase.Contains("Default")) mode=0;
4428 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4430 Int_t shieldindex=0;
4431 if (shield.Contains("SDD")) { // SDDouter
4432 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4434 } else if (shield.Contains("SPD")) { // SPDouter
4435 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4438 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4442 // do nothing if we are already beyond the shield
4443 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4444 if(dir<0 && rTrack > rToGo) return 1; // going outward
4445 if(dir>0 && rTrack < rToGo) return 1; // going inward
4449 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4451 Int_t index=2*fCurrentEsdTrack+shieldindex;
4453 Double_t xOverX0,x0,lengthTimesMeanDensity;
4458 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4459 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4460 lengthTimesMeanDensity = xOverX0*x0;
4461 lengthTimesMeanDensity *= dir;
4462 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4465 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4466 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4469 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4470 xOverX0 = fxOverX0Shield[shieldindex];
4471 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4472 lengthTimesMeanDensity *= dir;
4473 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4476 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4477 if(fxOverX0ShieldTrks[index]<0) {
4478 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4479 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4480 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4481 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4482 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4485 xOverX0 = fxOverX0ShieldTrks[index];
4486 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4487 lengthTimesMeanDensity *= dir;
4488 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4494 //------------------------------------------------------------------------
4495 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4497 Double_t oldGlobXYZ[3],
4498 TString direction) {
4499 //-------------------------------------------------------------------
4500 // Propagate beyond layer and correct for material
4501 // (material budget in different ways according to fUseTGeo value)
4502 // Add time if going outward (PropagateTo or PropagateToTGeo)
4503 //-------------------------------------------------------------------
4505 // Define budget mode:
4506 // 0: material from AliITSRecoParam (hard coded)
4507 // 1: material from TGeo in stepsof X cm (on the fly)
4508 // X = AliITSRecoParam::GetStepSizeTGeo()
4509 // 2: material from lut
4510 // 3: material from TGeo in one step (same for all hypotheses)
4523 if(fTrackingPhase.Contains("Clusters2Tracks"))
4524 { mode=3; } else { mode=1; }
4527 if(fTrackingPhase.Contains("Clusters2Tracks"))
4528 { mode=3; } else { mode=2; }
4534 if(fTrackingPhase.Contains("Default")) mode=0;
4536 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4538 Double_t r=fgLayers[layerindex].GetR();
4539 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4541 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4543 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4545 Int_t index=6*fCurrentEsdTrack+layerindex;
4548 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4551 // back before material (no correction)
4553 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4554 if (!t->GetLocalXat(rOld,xOld)) return 0;
4555 if (!t->Propagate(xOld)) return 0;
4559 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4560 lengthTimesMeanDensity = xOverX0*x0;
4561 lengthTimesMeanDensity *= dir;
4562 // Bring the track beyond the material
4563 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4566 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4567 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4570 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4571 xOverX0 = fxOverX0Layer[layerindex];
4572 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4573 lengthTimesMeanDensity *= dir;
4574 // Bring the track beyond the material
4575 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4578 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4579 if(fxOverX0LayerTrks[index]<0) {
4580 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4581 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4582 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4583 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4584 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4585 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4588 xOverX0 = fxOverX0LayerTrks[index];
4589 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4590 lengthTimesMeanDensity *= dir;
4591 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4598 //------------------------------------------------------------------------
4599 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4600 //-----------------------------------------------------------------
4601 // Initialize LUT for storing material for each prolonged track
4602 //-----------------------------------------------------------------
4603 fxOverX0PipeTrks = new Float_t[ntracks];
4604 fxTimesRhoPipeTrks = new Float_t[ntracks];
4605 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4606 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4607 fxOverX0LayerTrks = new Float_t[ntracks*6];
4608 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4610 for(Int_t i=0; i<ntracks; i++) {
4611 fxOverX0PipeTrks[i] = -1.;
4612 fxTimesRhoPipeTrks[i] = -1.;
4614 for(Int_t j=0; j<ntracks*2; j++) {
4615 fxOverX0ShieldTrks[j] = -1.;
4616 fxTimesRhoShieldTrks[j] = -1.;
4618 for(Int_t k=0; k<ntracks*6; k++) {
4619 fxOverX0LayerTrks[k] = -1.;
4620 fxTimesRhoLayerTrks[k] = -1.;
4627 //------------------------------------------------------------------------
4628 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4629 //-----------------------------------------------------------------
4630 // Delete LUT for storing material for each prolonged track
4631 //-----------------------------------------------------------------
4632 if(fxOverX0PipeTrks) {
4633 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4635 if(fxOverX0ShieldTrks) {
4636 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4639 if(fxOverX0LayerTrks) {
4640 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4642 if(fxTimesRhoPipeTrks) {
4643 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4645 if(fxTimesRhoShieldTrks) {
4646 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4648 if(fxTimesRhoLayerTrks) {
4649 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4653 //------------------------------------------------------------------------
4654 void AliITStrackerMI::SetForceSkippingOfLayer() {
4655 //-----------------------------------------------------------------
4656 // Check if we are forced to skip layers
4657 // either we set to skip them in RecoParam
4658 // or they were off during data-taking
4659 //-----------------------------------------------------------------
4661 const AliEventInfo *eventInfo = GetEventInfo();
4663 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4664 fForceSkippingOfLayer[l] = 0;
4666 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4670 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4671 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4673 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4674 } else if(l==2 || l==3) {
4675 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4677 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4683 //------------------------------------------------------------------------
4684 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4685 Int_t ilayer,Int_t idet) const {
4686 //-----------------------------------------------------------------
4687 // This method is used to decide whether to allow a prolongation
4688 // without clusters, because we want to skip the layer.
4689 // In this case the return value is > 0:
4690 // return 1: the user requested to skip a layer
4691 // return 2: track outside z acceptance
4692 //-----------------------------------------------------------------
4694 if (ForceSkippingOfLayer(ilayer)) return 1;
4696 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4698 if (idet<0 && // out in z
4699 ilayer>innerLayCanSkip &&
4700 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4701 // check if track will cross SPD outer layer
4702 Double_t phiAtSPD2,zAtSPD2;
4703 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4704 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4706 return 2; // always allow skipping, changed on 05.11.2009
4711 //------------------------------------------------------------------------
4712 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4713 Int_t ilayer,Int_t idet,
4714 Double_t dz,Double_t dy,
4715 Bool_t noClusters) const {
4716 //-----------------------------------------------------------------
4717 // This method is used to decide whether to allow a prolongation
4718 // without clusters, because there is a dead zone in the road.
4719 // In this case the return value is > 0:
4720 // return 1: dead zone at z=0,+-7cm in SPD
4721 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4722 // return 2: all road is "bad" (dead or noisy) from the OCDB
4723 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4724 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4725 //-----------------------------------------------------------------
4727 // check dead zones at z=0,+-7cm in the SPD
4728 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4729 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4730 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4731 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4732 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4733 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4734 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4735 for (Int_t i=0; i<3; i++)
4736 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4737 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4738 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4742 // check bad zones from OCDB
4743 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4745 if (idet<0) return 0;
4747 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4750 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4751 if (ilayer==0 || ilayer==1) { // ---------- SPD
4753 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4755 detSizeFactorX *= 2.;
4756 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4759 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4760 if (detType==2) segm->SetLayer(ilayer+1);
4761 Float_t detSizeX = detSizeFactorX*segm->Dx();
4762 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4764 // check if the road overlaps with bad chips
4766 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4767 Float_t zlocmin = zloc-dz;
4768 Float_t zlocmax = zloc+dz;
4769 Float_t xlocmin = xloc-dy;
4770 Float_t xlocmax = xloc+dy;
4771 Int_t chipsInRoad[100];
4773 // check if road goes out of detector
4774 Bool_t touchNeighbourDet=kFALSE;
4775 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4776 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4777 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4778 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4779 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));
4781 // check if this detector is bad
4783 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4784 if(!touchNeighbourDet) {
4785 return 2; // all detectors in road are bad
4787 return 3; // at least one is bad
4791 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4792 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4793 if (!nChipsInRoad) return 0;
4795 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4796 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4797 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4798 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4799 if (det.IsChipBad(chipsInRoad[iCh])) {
4807 if(!touchNeighbourDet) {
4808 AliDebug(2,"all bad in road");
4809 return 2; // all chips in road are bad
4811 return 3; // at least a bad chip in road
4816 AliDebug(2,"at least a bad in road");
4817 return 3; // at least a bad chip in road
4821 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4822 || !noClusters) return 0;
4824 // There are no clusters in road: check if there is at least
4825 // a bad SPD pixel or SDD anode or SSD strips on both sides
4827 Int_t idetInITS=idet;
4828 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4830 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4831 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4834 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4838 //------------------------------------------------------------------------
4839 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4840 const AliITStrackMI *track,
4841 Float_t &xloc,Float_t &zloc) const {
4842 //-----------------------------------------------------------------
4843 // Gives position of track in local module ref. frame
4844 //-----------------------------------------------------------------
4849 if(idet<0) return kTRUE; // track out of z acceptance of layer
4851 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4853 Int_t lad = Int_t(idet/ndet) + 1;
4855 Int_t det = idet - (lad-1)*ndet + 1;
4857 Double_t xyzGlob[3],xyzLoc[3];
4859 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4860 // take into account the misalignment: xyz at real detector plane
4861 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4863 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4865 xloc = (Float_t)xyzLoc[0];
4866 zloc = (Float_t)xyzLoc[2];
4870 //------------------------------------------------------------------------
4871 //------------------------------------------------------------------------
4872 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4874 // Method to be optimized further:
4875 // Aim: decide whether a track can be used for PlaneEff evaluation
4876 // the decision is taken based on the track quality at the layer under study
4877 // no information on the clusters on this layer has to be used
4878 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4879 // the cut is done on number of sigmas from the boundaries
4881 // Input: Actual track, layer [0,5] under study
4883 // Return: kTRUE if this is a good track
4885 // it will apply a pre-selection to obtain good quality tracks.
4886 // Here also you will have the possibility to put a control on the
4887 // impact point of the track on the basic block, in order to exclude border regions
4888 // this will be done by calling a proper method of the AliITSPlaneEff class.
4890 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4891 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4893 Int_t index[AliITSgeomTGeo::kNLayers];
4895 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4897 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4898 index[k]=clusters[k];
4902 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4903 AliITSlayer &layer=fgLayers[ilayer];
4904 Double_t r=layer.GetR();
4905 AliITStrackMI tmp(*track);
4907 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4908 Int_t ncl_out=0; Int_t ncl_in=0;
4909 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4910 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4911 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4912 // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4913 if(index[lay]>=0)ncl_out++;
4915 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4916 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4917 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4918 if (index[lay]>=0) ncl_in++;
4920 Int_t ncl=ncl_out+ncl_in;
4921 Bool_t nextout = kFALSE;
4922 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4923 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4924 Bool_t nextin = kFALSE;
4925 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4926 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4927 // maximum number of missing clusters allowed in outermost layers
4928 if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4930 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4931 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4933 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4934 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4935 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4936 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4940 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4941 Int_t idet=layer.FindDetectorIndex(phi,z);
4942 if(idet<0) { AliInfo(Form("cannot find detector"));
4945 // here check if it has good Chi Square.
4947 //propagate to the intersection with the detector plane
4948 const AliITSdetector &det=layer.GetDetector(idet);
4949 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4953 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4954 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4955 if(key>fPlaneEff->Nblock()) return kFALSE;
4956 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4957 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4959 // DEFINITION OF SEARCH ROAD FOR accepting a track
4961 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4962 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4963 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4964 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4965 // done for RecPoints
4967 // exclude tracks at boundary between detectors
4968 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4969 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4970 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4971 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4972 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4973 if ( (locx-dx < blockXmn+boundaryWidth) ||
4974 (locx+dx > blockXmx-boundaryWidth) ||
4975 (locz-dz < blockZmn+boundaryWidth) ||
4976 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4979 //------------------------------------------------------------------------
4980 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4982 // This Method has to be optimized! For the time-being it uses the same criteria
4983 // as those used in the search of extra clusters for overlapping modules.
4985 // Method Purpose: estabilish whether a track has produced a recpoint or not
4986 // in the layer under study (For Plane efficiency)
4988 // inputs: AliITStrackMI* track (pointer to a usable track)
4990 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4991 // traversed by this very track. In details:
4992 // - if a cluster can be associated to the track then call
4993 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4994 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4997 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4998 AliITSlayer &layer=fgLayers[ilayer];
4999 Double_t r=layer.GetR();
5000 AliITStrackMI tmp(*track);
5004 if (!tmp.GetPhiZat(r,phi,z)) return;
5005 Int_t idet=layer.FindDetectorIndex(phi,z);
5007 if(idet<0) { AliInfo(Form("cannot find detector"));
5011 //propagate to the intersection with the detector plane
5012 const AliITSdetector &det=layer.GetDetector(idet);
5013 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5017 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5019 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5020 TMath::Sqrt(tmp.GetSigmaZ2() +
5021 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5022 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5023 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5024 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5025 TMath::Sqrt(tmp.GetSigmaY2() +
5026 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5027 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5028 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5030 // road in global (rphi,z) [i.e. in tracking ref. system]
5031 Double_t zmin = tmp.GetZ() - dz;
5032 Double_t zmax = tmp.GetZ() + dz;
5033 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5034 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5036 // select clusters in road
5037 layer.SelectClusters(zmin,zmax,ymin,ymax);
5039 // Define criteria for track-cluster association
5040 Double_t msz = tmp.GetSigmaZ2() +
5041 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5042 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5043 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5044 Double_t msy = tmp.GetSigmaY2() +
5045 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5046 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5047 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5048 if (tmp.GetConstrain()) {
5049 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5050 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5052 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5053 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5055 msz = 1./msz; // 1/RoadZ^2
5056 msy = 1./msy; // 1/RoadY^2
5059 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5061 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5062 //Double_t tolerance=0.2;
5063 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5064 idetc = cl->GetDetectorIndex();
5065 if(idet!=idetc) continue;
5066 //Int_t ilay = cl->GetLayer();
5068 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5069 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5071 Double_t chi2=tmp.GetPredictedChi2(cl);
5072 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5076 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5078 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5079 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5080 if(key>fPlaneEff->Nblock()) return;
5081 Bool_t found=kFALSE;
5084 while ((cl=layer.GetNextCluster(clidx))!=0) {
5085 idetc = cl->GetDetectorIndex();
5086 if(idet!=idetc) continue;
5087 // here real control to see whether the cluster can be associated to the track.
5088 // cluster not associated to track
5089 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5090 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5091 // calculate track-clusters chi2
5092 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5093 // in particular, the error associated to the cluster
5094 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5096 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5098 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5099 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5100 // track->SetExtraModule(ilayer,idetExtra);
5102 if(!fPlaneEff->UpDatePlaneEff(found,key))
5103 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5105 // this for FO efficiency studies (only for SPD) //
5106 UInt_t keyFO=999999;
5107 Bool_t foundFO=kFALSE;
5108 if(ilayer<2){ //ONLY SPD layers for FastOr studies
5109 TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5110 Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5111 if(!fSPDChipIntPlaneEff[key]){
5112 AliITSPlaneEffSPD spd;
5113 keyFO = spd.SwitchChipKeyNumbering(key);
5114 if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5115 keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5116 if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5117 AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5120 if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5121 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5127 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5128 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5129 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5130 Int_t cltype[2]={-999,-999};
5133 Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5137 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5138 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5141 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5142 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5143 cltype[0]=layer.GetCluster(ci)->GetNy();
5144 cltype[1]=layer.GetCluster(ci)->GetNz();
5146 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5147 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5148 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5149 // It is computed properly by calling the method
5150 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5152 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5153 //if (tmp.PropagateTo(x,0.,0.)) {
5154 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5155 AliCluster c(*layer.GetCluster(ci));
5156 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5157 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5158 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5159 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5160 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5163 // Compute the angles between the track and the module
5164 // compute the angle "in phi direction", i.e. the angle in the transverse plane
5165 // between the normal to the module and the projection (in the transverse plane) of the
5167 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5168 Float_t tgl = tmp.GetTgl();
5169 Float_t phitr = tmp.GetSnp();
5170 phitr = TMath::ASin(phitr);
5171 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5173 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5174 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5176 alpha = tmp.GetAlpha();
5177 Double_t phiglob = alpha+phitr;
5179 p[0] = TMath::Cos(phiglob);
5180 p[1] = TMath::Sin(phiglob);
5182 TVector3 pvec(p[0],p[1],p[2]);
5183 TVector3 normvec(rot[1],rot[4],rot[7]);
5184 Double_t angle = pvec.Angle(normvec);
5186 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5187 angle *= 180./TMath::Pi();
5190 TVector3 pt(p[0],p[1],0);
5191 TVector3 normt(rot[1],rot[4],0);
5192 Double_t anglet = pt.Angle(normt);
5194 Double_t phiPt = TMath::ATan2(p[1],p[0]);
5195 if(phiPt<0)phiPt+=2.*TMath::Pi();
5196 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5197 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5198 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5199 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
5200 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5201 anglet *= 180./TMath::Pi();
5203 AngleModTrack[2]=(Float_t) angle;
5204 AngleModTrack[0]=(Float_t) anglet;
5205 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5206 AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5207 AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5208 AngleModTrack[1]*=180./TMath::Pi(); // in degree
5210 fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
5212 // For FO efficiency studies of SPD
5213 if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,AngleModTrack);
5215 if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5219 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5221 // find the MC cluster for the label
5222 return fgLayers[lr].FindClusterForLabel(label,store);
5226 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5228 // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5229 strncpy(patt,"......",6);
5231 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5234 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5235 Int_t cindex = track->GetClusterIndex(i);
5236 Int_t l=(cindex & 0xf0000000) >> 28;
5237 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5239 for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5240 patt[l] = isWrong ? 'f':'c';
5246 //------------------------------------------------------------------------
5247 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5249 //--------------------------------------------------------------------
5252 for (int ic=0;ic<fN;ic++) {
5253 const AliITSRecPoint *cl = GetCluster(ic);
5254 for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5256 if (store) store[nfound] = ic;