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<24;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 //------------------------------------------------------------------
3404 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3406 array = new TObjArray(10);
3407 fTrackHypothesys.AddAtAndExpand(array,esdindex);
3409 array->AddLast(track);
3411 //------------------------------------------------------------------------
3412 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3414 //-------------------------------------------------------------------
3415 // compress array of track hypothesys
3416 // keep only maxsize best hypothesys
3417 //-------------------------------------------------------------------
3418 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3419 if (! (fTrackHypothesys.At(esdindex)) ) return;
3420 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3421 Int_t entries = array->GetEntriesFast();
3423 //- find preliminary besttrack as a reference
3424 Float_t minchi2=10000;
3426 AliITStrackMI * besttrack=0;
3428 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3429 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3430 if (!track) continue;
3431 Float_t chi2 = NormalizedChi2(track,0);
3433 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3434 track->SetLabel(tpcLabel);
3436 track->SetFakeRatio(1.);
3437 CookLabel(track,0.); //For comparison only
3439 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3440 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3441 if (track->GetNumberOfClusters()<maxn) continue;
3442 maxn = track->GetNumberOfClusters();
3443 // if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3450 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3451 delete array->RemoveAt(itrack);
3455 if (!besttrack) return;
3458 //take errors of best track as a reference
3459 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3460 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3461 for (Int_t j=0;j<6;j++) {
3462 if (besttrack->GetClIndex(j)>=0){
3463 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3464 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3465 ny[j] = besttrack->GetNy(j);
3466 nz[j] = besttrack->GetNz(j);
3470 // calculate normalized chi2
3472 Float_t * chi2 = new Float_t[entries];
3473 Int_t * index = new Int_t[entries];
3474 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3475 for (Int_t itrack=0;itrack<entries;itrack++){
3476 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3478 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3479 double chi2t = GetNormalizedChi2(track, mode);
3480 track->SetChi2MIP(0,chi2t);
3481 if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3482 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3483 chi2[itrack] = chi2t;
3486 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3487 delete array->RemoveAt(itrack);
3493 TMath::Sort(entries,chi2,index,kFALSE);
3494 besttrack = (AliITStrackMI*)array->At(index[0]);
3495 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3496 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3497 for (Int_t j=0;j<6;j++){
3498 if (besttrack->GetClIndex(j)>=0){
3499 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3500 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3501 ny[j] = besttrack->GetNy(j);
3502 nz[j] = besttrack->GetNz(j);
3507 // calculate one more time with updated normalized errors
3508 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3509 for (Int_t itrack=0;itrack<entries;itrack++){
3510 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3512 double chi2t = GetNormalizedChi2(track, mode);
3513 track->SetChi2MIP(0,chi2t);
3514 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3515 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3516 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3517 chi2[itrack] = chi2t; //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3520 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3521 delete array->RemoveAt(itrack);
3526 entries = array->GetEntriesFast();
3530 TObjArray * newarray = new TObjArray();
3531 TMath::Sort(entries,chi2,index,kFALSE);
3532 besttrack = (AliITStrackMI*)array->At(index[0]);
3534 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3536 for (Int_t j=0;j<6;j++){
3537 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3538 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3539 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3540 ny[j] = besttrack->GetNy(j);
3541 nz[j] = besttrack->GetNz(j);
3544 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3545 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3546 Float_t minn = besttrack->GetNumberOfClusters()-3;
3548 for (Int_t i=0;i<entries;i++){
3549 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3550 if (!track) continue;
3551 if (accepted>maxcut) break;
3552 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3553 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3554 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3555 delete array->RemoveAt(index[i]);
3559 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3560 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3561 if (!shortbest) accepted++;
3563 newarray->AddLast(array->RemoveAt(index[i]));
3564 for (Int_t j=0;j<6;j++){
3566 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3567 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3568 ny[j] = track->GetNy(j);
3569 nz[j] = track->GetNz(j);
3574 delete array->RemoveAt(index[i]);
3578 delete fTrackHypothesys.RemoveAt(esdindex);
3579 fTrackHypothesys.AddAt(newarray,esdindex);
3583 delete fTrackHypothesys.RemoveAt(esdindex);
3589 //------------------------------------------------------------------------
3590 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3592 //-------------------------------------------------------------
3593 // try to find best hypothesy
3594 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3595 // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3596 //-------------------------------------------------------------
3597 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3598 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3599 if (!array) return 0;
3600 Int_t entries = array->GetEntriesFast();
3601 if (!entries) return 0;
3602 Float_t minchi2 = 100000;
3603 AliITStrackMI * besttrack=0;
3605 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3606 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3607 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3608 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3610 for (Int_t i=0;i<entries;i++){
3611 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3612 if (!track) continue;
3613 Float_t sigmarfi,sigmaz;
3614 GetDCASigma(track,sigmarfi,sigmaz);
3615 track->SetDnorm(0,sigmarfi);
3616 track->SetDnorm(1,sigmaz);
3618 track->SetChi2MIP(1,1000000);
3619 track->SetChi2MIP(2,1000000);
3620 track->SetChi2MIP(3,1000000);
3623 backtrack = new(backtrack) AliITStrackMI(*track);
3624 if (track->GetConstrain()) {
3625 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3626 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3627 if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3628 else {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3630 backtrack->ResetCovariance(10.);
3632 backtrack->ResetCovariance(10.);
3634 backtrack->ResetClusters();
3636 Double_t x = original->GetX();
3637 if (!RefitAt(x,backtrack,track)) continue;
3639 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3640 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3641 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3642 track->SetChi22(GetMatchingChi2(backtrack,original));
3644 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3645 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3646 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3649 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3651 //forward track - without constraint
3652 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3653 forwardtrack->ResetClusters();
3655 if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue; // w/o fwd track MIP03 is meaningless
3656 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3657 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3658 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3660 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3661 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3662 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3663 forwardtrack->SetD(0,track->GetD(0));
3664 forwardtrack->SetD(1,track->GetD(1));
3667 AliITSRecPoint* clist[6];
3668 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3669 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3672 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3673 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3674 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3675 track->SetChi2MIP(3,1000);
3678 Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed(); //RS
3679 if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3680 else chi2 += track->GetNUsed();
3682 for (Int_t ichi=0;ichi<5;ichi++){
3683 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3685 if (chi2 < minchi2){
3686 //besttrack = new AliITStrackMI(*forwardtrack);
3688 besttrack->SetLabel(track->GetLabel());
3689 besttrack->SetFakeRatio(track->GetFakeRatio());
3691 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3692 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3693 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3697 delete forwardtrack;
3699 if (!besttrack) return 0;
3702 for (Int_t i=0;i<entries;i++){
3703 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3705 if (!track) continue;
3707 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3708 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3709 // RS: don't apply this cut when fSelectBestMIP03 is on
3710 || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3712 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3713 delete array->RemoveAt(i);
3723 SortTrackHypothesys(esdindex,checkmax,1);
3725 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3726 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3727 besttrack = (AliITStrackMI*)array->At(0);
3728 if (!besttrack) return 0;
3729 besttrack->SetChi2MIP(8,0);
3730 fBestTrackIndex[esdindex]=0;
3731 entries = array->GetEntriesFast();
3732 AliITStrackMI *longtrack =0;
3734 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3735 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3736 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3737 if (!track->GetConstrain()) continue;
3738 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3739 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3740 if (track->GetChi2MIP(0)>4.) continue;
3741 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3744 //if (longtrack) besttrack=longtrack;
3746 // RS do shared cluster analysis here only if the new sharing analysis is not requested
3747 //RRR if (fFlagFakes) return besttrack;
3750 AliITSRecPoint * clist[6];
3751 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3752 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3753 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3754 RegisterClusterTracks(besttrack,esdindex);
3761 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3762 if (sharedtrack>=0){
3764 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3766 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3772 if (shared>2.5) return 0;
3773 if (shared>1.0) return besttrack;
3775 // Don't sign clusters if not gold track
3777 if (!besttrack->IsGoldPrimary()) return besttrack;
3778 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3780 if (fConstraint[fPass]){
3784 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3785 for (Int_t i=0;i<6;i++){
3786 Int_t index = besttrack->GetClIndex(i);
3787 if (index<0) continue;
3788 Int_t ilayer = (index & 0xf0000000) >> 28;
3789 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3790 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3792 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3793 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3794 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3795 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3796 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3797 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3799 Bool_t cansign = kTRUE;
3800 for (Int_t itrack=0;itrack<entries; itrack++){
3801 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3802 if (!track) continue;
3803 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3804 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3810 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3811 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3812 if (!c->IsUsed()) c->Use();
3818 //------------------------------------------------------------------------
3819 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3822 // get "best" hypothesys
3825 Int_t nentries = itsTracks.GetEntriesFast();
3826 for (Int_t i=0;i<nentries;i++){
3827 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3828 if (!track) continue;
3829 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3830 if (!array) continue;
3831 if (array->GetEntriesFast()<=0) continue;
3833 AliITStrackMI* longtrack=0;
3835 Float_t maxchi2=1000;
3836 for (Int_t j=0;j<array->GetEntriesFast();j++){
3837 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3838 if (!trackHyp) continue;
3839 if (trackHyp->GetGoldV0()) {
3840 longtrack = trackHyp; //gold V0 track taken
3843 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3844 Float_t chi2 = trackHyp->GetChi2MIP(0);
3845 if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3846 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3848 if (fAfterV0){ // ??? RS
3849 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3851 if (chi2 > maxchi2) continue;
3852 minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3853 if (fSelectBestMIP03) minn++; // allow next to longest to win
3860 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3861 if (!longtrack) {longtrack = besttrack;}
3862 else besttrack= longtrack;
3866 AliITSRecPoint * clist[6];
3867 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3869 track->SetNUsed(shared);
3870 track->SetNSkipped(besttrack->GetNSkipped());
3871 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3873 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3877 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3878 //if (sharedtrack==-1) sharedtrack=0;
3879 if (sharedtrack>=0) {
3880 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3883 if (besttrack&&fAfterV0) {
3884 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3885 track->SetWinner(besttrack);
3888 if (fConstraint[fPass]) {
3889 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3890 track->SetWinner(besttrack);
3892 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3893 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3894 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3901 //------------------------------------------------------------------------
3902 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3905 // RS: flag those tracks which are suxpected to have fake clusters
3907 const double kThreshPt = 0.5;
3908 AliRefArray *refArr[6];
3910 for (int i=0;i<6;i++) {
3911 int ncl = fgLayers[i].GetNumberOfClusters();
3912 refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3914 Int_t nentries = itsTracks.GetEntriesFast();
3916 // fill cluster->track associations
3917 for (Int_t itr=0;itr<nentries;itr++){
3918 AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3919 if (!track) continue;
3920 AliITStrackMI* trackITS = track->GetWinner();
3921 if (!trackITS) continue;
3922 for (int il=trackITS->GetNumberOfClusters();il--;) {
3923 int idx = trackITS->GetClusterIndex(il);
3924 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3925 // if (c>fgLayers[l].GetNumberOfClusters()) continue;
3926 refArr[l]->AddReference(c, itr);
3930 const UInt_t kMaxRef = 100;
3931 UInt_t crefs[kMaxRef];
3933 // process tracks with shared clusters
3934 for (int itr=0;itr<nentries;itr++){
3935 AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3936 AliITStrackMI* trackH0 = track0->GetWinner();
3937 if (!trackH0) continue;
3938 AliESDtrack* esd0 = track0->GetESDtrack();
3940 for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3941 int idx = trackH0->GetClusterIndex(il);
3942 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3943 ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3944 if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3945 esd0->SetITSSharedFlag(l);
3946 for (int ir=ncrefs;ir--;) {
3947 if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3948 AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3949 AliITStrackMI* trackH1 = track1->GetWinner();
3950 AliESDtrack* esd1 = track1->GetESDtrack();
3951 esd1->SetITSSharedFlag(l);
3953 double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3954 if (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3955 else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3957 // select the one with smallest chi2's product
3958 res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3959 res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3961 if (res<0) esd1->SetITSFakeFlag(); // esd0 is winner
3962 else esd0->SetITSFakeFlag(); // esd1 is winner
3969 for (int i=6;i--;) delete refArr[i];
3972 //------------------------------------------------------------------------
3973 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3974 //--------------------------------------------------------------------
3975 //This function "cooks" a track label. If label<0, this track is fake.
3976 //--------------------------------------------------------------------
3979 if (track->GetESDtrack()){
3980 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3981 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3982 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3984 track->SetChi2MIP(9,0);
3986 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3987 Int_t cindex = track->GetClusterIndex(i);
3988 Int_t l=(cindex & 0xf0000000) >> 28;
3989 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3991 for (Int_t ind=0;ind<3;ind++){
3992 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3993 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3995 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3998 Int_t nclusters = track->GetNumberOfClusters();
3999 if (nclusters > 0) //PH Some tracks don't have any cluster
4000 track->SetFakeRatio(double(nwrong)/double(nclusters));
4001 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4002 track->SetLabel(-tpcLabel);
4004 track->SetLabel(tpcLabel);
4006 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4009 //------------------------------------------------------------------------
4010 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4012 // Fill the dE/dx in this track
4014 track->SetChi2MIP(9,0);
4015 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4016 Int_t cindex = track->GetClusterIndex(i);
4017 Int_t l=(cindex & 0xf0000000) >> 28;
4018 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4019 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4021 for (Int_t ind=0;ind<3;ind++){
4022 if (cl->GetLabel(ind)==lab) isWrong=0;
4024 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4028 track->CookdEdx(low,up);
4030 //------------------------------------------------------------------------
4031 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4033 // Create some arrays
4035 if (fCoefficients) delete []fCoefficients;
4036 fCoefficients = new Float_t[ntracks*48];
4037 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4039 //------------------------------------------------------------------------
4040 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4043 // Compute predicted chi2
4045 // Take into account the mis-alignment (bring track to cluster plane)
4046 Double_t xTrOrig=track->GetX();
4047 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4048 Float_t erry,errz,covyz;
4049 Float_t theta = track->GetTgl();
4050 Float_t phi = track->GetSnp();
4051 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4052 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4053 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()));
4054 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()));
4055 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4056 // Bring the track back to detector plane in ideal geometry
4057 // [mis-alignment will be accounted for in UpdateMI()]
4058 if (!track->Propagate(xTrOrig)) return 1000.;
4060 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4061 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4063 chi2+=0.5*TMath::Min(delta/2,2.);
4064 chi2+=2.*cluster->GetDeltaProbability();
4067 track->SetNy(layer,ny);
4068 track->SetNz(layer,nz);
4069 track->SetSigmaY(layer,erry);
4070 track->SetSigmaZ(layer, errz);
4071 track->SetSigmaYZ(layer,covyz);
4072 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4073 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4077 //------------------------------------------------------------------------
4078 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4083 Int_t layer = (index & 0xf0000000) >> 28;
4084 track->SetClIndex(layer, index);
4085 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4086 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4087 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4088 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4092 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
4095 // Take into account the mis-alignment (bring track to cluster plane)
4096 Double_t xTrOrig=track->GetX();
4097 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4098 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4099 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4100 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
4102 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4105 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4106 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4107 c.SetSigmaYZ(track->GetSigmaYZ(layer));
4110 Int_t updated = track->UpdateMI(&c,chi2,index);
4111 // Bring the track back to detector plane in ideal geometry
4112 if (!track->Propagate(xTrOrig)) return 0;
4114 if(!updated) AliDebug(2,"update failed");
4118 //------------------------------------------------------------------------
4119 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4122 //DCA sigmas parameterization
4123 //to be paramterized using external parameters in future
4126 Double_t curv=track->GetC();
4127 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4128 sigmaz = 0.0110+4.37*TMath::Abs(curv);
4130 //------------------------------------------------------------------------
4131 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4134 // Clusters from delta electrons?
4136 Int_t entries = clusterArray->GetEntriesFast();
4137 if (entries<4) return;
4138 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4139 Int_t layer = cluster->GetLayer();
4140 if (layer>1) return;
4142 Int_t ncandidates=0;
4143 Float_t r = (layer>0)? 7:4;
4145 for (Int_t i=0;i<entries;i++){
4146 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4147 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4148 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4149 index[ncandidates] = i; //candidate to belong to delta electron track
4151 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4152 cl0->SetDeltaProbability(1);
4158 for (Int_t i=0;i<ncandidates;i++){
4159 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4160 if (cl0->GetDeltaProbability()>0.8) continue;
4163 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4164 sumy=sumz=sumy2=sumyz=sumw=0.0;
4165 for (Int_t j=0;j<ncandidates;j++){
4167 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4169 Float_t dz = cl0->GetZ()-cl1->GetZ();
4170 Float_t dy = cl0->GetY()-cl1->GetY();
4171 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4172 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4173 y[ncl] = cl1->GetY();
4174 z[ncl] = cl1->GetZ();
4175 sumy+= y[ncl]*weight;
4176 sumz+= z[ncl]*weight;
4177 sumy2+=y[ncl]*y[ncl]*weight;
4178 sumyz+=y[ncl]*z[ncl]*weight;
4183 if (ncl<4) continue;
4184 Float_t det = sumw*sumy2 - sumy*sumy;
4186 if (TMath::Abs(det)>0.01){
4187 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
4188 Float_t k = (sumyz*sumw - sumy*sumz)/det;
4189 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4192 Float_t z0 = sumyz/sumy;
4193 delta = TMath::Abs(cl0->GetZ()-z0);
4196 cl0->SetDeltaProbability(1-20.*delta);
4200 //------------------------------------------------------------------------
4201 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4206 track->UpdateESDtrack(flags);
4207 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4208 if (oldtrack) delete oldtrack;
4209 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4210 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4211 // printf("Problem\n");
4214 //------------------------------------------------------------------------
4215 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4217 // Get nearest upper layer close to the point xr.
4218 // rough approximation
4220 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4221 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4223 for (Int_t i=0;i<6;i++){
4224 if (radius<kRadiuses[i]){
4231 //------------------------------------------------------------------------
4232 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4233 //--------------------------------------------------------------------
4234 // Fill a look-up table with mean material
4235 //--------------------------------------------------------------------
4239 Double_t point1[3],point2[3];
4240 Double_t phi,cosphi,sinphi,z;
4241 // 0-5 layers, 6 pipe, 7-8 shields
4242 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4243 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4245 Int_t ifirst=0,ilast=0;
4246 if(material.Contains("Pipe")) {
4248 } else if(material.Contains("Shields")) {
4250 } else if(material.Contains("Layers")) {
4253 Error("BuildMaterialLUT","Wrong layer name\n");
4256 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4257 Double_t param[5]={0.,0.,0.,0.,0.};
4258 for (Int_t i=0; i<n; i++) {
4259 phi = 2.*TMath::Pi()*gRandom->Rndm();
4260 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4261 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4262 point1[0] = rmin[imat]*cosphi;
4263 point1[1] = rmin[imat]*sinphi;
4265 point2[0] = rmax[imat]*cosphi;
4266 point2[1] = rmax[imat]*sinphi;
4268 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4269 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4271 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4273 fxOverX0Layer[imat] = param[1];
4274 fxTimesRhoLayer[imat] = param[0]*param[4];
4275 } else if(imat==6) {
4276 fxOverX0Pipe = param[1];
4277 fxTimesRhoPipe = param[0]*param[4];
4278 } else if(imat==7) {
4279 fxOverX0Shield[0] = param[1];
4280 fxTimesRhoShield[0] = param[0]*param[4];
4281 } else if(imat==8) {
4282 fxOverX0Shield[1] = param[1];
4283 fxTimesRhoShield[1] = param[0]*param[4];
4287 printf("%s\n",material.Data());
4288 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4289 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4290 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4291 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4292 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4293 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4294 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4295 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4296 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4300 //------------------------------------------------------------------------
4301 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4302 TString direction) {
4303 //-------------------------------------------------------------------
4304 // Propagate beyond beam pipe and correct for material
4305 // (material budget in different ways according to fUseTGeo value)
4306 // Add time if going outward (PropagateTo or PropagateToTGeo)
4307 //-------------------------------------------------------------------
4309 // Define budget mode:
4310 // 0: material from AliITSRecoParam (hard coded)
4311 // 1: material from TGeo in one step (on the fly)
4312 // 2: material from lut
4313 // 3: material from TGeo in one step (same for all hypotheses)
4326 if(fTrackingPhase.Contains("Clusters2Tracks"))
4327 { mode=3; } else { mode=1; }
4330 if(fTrackingPhase.Contains("Clusters2Tracks"))
4331 { mode=3; } else { mode=2; }
4337 if(fTrackingPhase.Contains("Default")) mode=0;
4339 Int_t index=fCurrentEsdTrack;
4341 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4342 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4344 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4346 Double_t xOverX0,x0,lengthTimesMeanDensity;
4350 xOverX0 = AliITSRecoParam::GetdPipe();
4351 x0 = AliITSRecoParam::GetX0Be();
4352 lengthTimesMeanDensity = xOverX0*x0;
4353 lengthTimesMeanDensity *= dir;
4354 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4357 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4360 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4361 xOverX0 = fxOverX0Pipe;
4362 lengthTimesMeanDensity = fxTimesRhoPipe;
4363 lengthTimesMeanDensity *= dir;
4364 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4367 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4368 if(fxOverX0PipeTrks[index]<0) {
4369 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4370 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4371 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4372 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4373 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4376 xOverX0 = fxOverX0PipeTrks[index];
4377 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4378 lengthTimesMeanDensity *= dir;
4379 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4385 //------------------------------------------------------------------------
4386 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4388 TString direction) {
4389 //-------------------------------------------------------------------
4390 // Propagate beyond SPD or SDD shield and correct for material
4391 // (material budget in different ways according to fUseTGeo value)
4392 // Add time if going outward (PropagateTo or PropagateToTGeo)
4393 //-------------------------------------------------------------------
4395 // Define budget mode:
4396 // 0: material from AliITSRecoParam (hard coded)
4397 // 1: material from TGeo in steps of X cm (on the fly)
4398 // X = AliITSRecoParam::GetStepSizeTGeo()
4399 // 2: material from lut
4400 // 3: material from TGeo in one step (same for all hypotheses)
4413 if(fTrackingPhase.Contains("Clusters2Tracks"))
4414 { mode=3; } else { mode=1; }
4417 if(fTrackingPhase.Contains("Clusters2Tracks"))
4418 { mode=3; } else { mode=2; }
4424 if(fTrackingPhase.Contains("Default")) mode=0;
4426 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4428 Int_t shieldindex=0;
4429 if (shield.Contains("SDD")) { // SDDouter
4430 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4432 } else if (shield.Contains("SPD")) { // SPDouter
4433 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4436 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4440 // do nothing if we are already beyond the shield
4441 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4442 if(dir<0 && rTrack > rToGo) return 1; // going outward
4443 if(dir>0 && rTrack < rToGo) return 1; // going inward
4447 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4449 Int_t index=2*fCurrentEsdTrack+shieldindex;
4451 Double_t xOverX0,x0,lengthTimesMeanDensity;
4456 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4457 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4458 lengthTimesMeanDensity = xOverX0*x0;
4459 lengthTimesMeanDensity *= dir;
4460 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4463 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4464 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4467 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4468 xOverX0 = fxOverX0Shield[shieldindex];
4469 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4470 lengthTimesMeanDensity *= dir;
4471 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4474 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4475 if(fxOverX0ShieldTrks[index]<0) {
4476 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4477 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4478 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4479 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4480 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4483 xOverX0 = fxOverX0ShieldTrks[index];
4484 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4485 lengthTimesMeanDensity *= dir;
4486 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4492 //------------------------------------------------------------------------
4493 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4495 Double_t oldGlobXYZ[3],
4496 TString direction) {
4497 //-------------------------------------------------------------------
4498 // Propagate beyond layer and correct for material
4499 // (material budget in different ways according to fUseTGeo value)
4500 // Add time if going outward (PropagateTo or PropagateToTGeo)
4501 //-------------------------------------------------------------------
4503 // Define budget mode:
4504 // 0: material from AliITSRecoParam (hard coded)
4505 // 1: material from TGeo in stepsof X cm (on the fly)
4506 // X = AliITSRecoParam::GetStepSizeTGeo()
4507 // 2: material from lut
4508 // 3: material from TGeo in one step (same for all hypotheses)
4521 if(fTrackingPhase.Contains("Clusters2Tracks"))
4522 { mode=3; } else { mode=1; }
4525 if(fTrackingPhase.Contains("Clusters2Tracks"))
4526 { mode=3; } else { mode=2; }
4532 if(fTrackingPhase.Contains("Default")) mode=0;
4534 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4536 Double_t r=fgLayers[layerindex].GetR();
4537 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4539 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4541 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4543 Int_t index=6*fCurrentEsdTrack+layerindex;
4546 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4549 // back before material (no correction)
4551 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4552 if (!t->GetLocalXat(rOld,xOld)) return 0;
4553 if (!t->Propagate(xOld)) return 0;
4557 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4558 lengthTimesMeanDensity = xOverX0*x0;
4559 lengthTimesMeanDensity *= dir;
4560 // Bring the track beyond the material
4561 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4564 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4565 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4568 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4569 xOverX0 = fxOverX0Layer[layerindex];
4570 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4571 lengthTimesMeanDensity *= dir;
4572 // Bring the track beyond the material
4573 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4576 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4577 if(fxOverX0LayerTrks[index]<0) {
4578 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4579 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4580 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4581 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4582 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4583 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4586 xOverX0 = fxOverX0LayerTrks[index];
4587 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4588 lengthTimesMeanDensity *= dir;
4589 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4596 //------------------------------------------------------------------------
4597 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4598 //-----------------------------------------------------------------
4599 // Initialize LUT for storing material for each prolonged track
4600 //-----------------------------------------------------------------
4601 fxOverX0PipeTrks = new Float_t[ntracks];
4602 fxTimesRhoPipeTrks = new Float_t[ntracks];
4603 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4604 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4605 fxOverX0LayerTrks = new Float_t[ntracks*6];
4606 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4608 for(Int_t i=0; i<ntracks; i++) {
4609 fxOverX0PipeTrks[i] = -1.;
4610 fxTimesRhoPipeTrks[i] = -1.;
4612 for(Int_t j=0; j<ntracks*2; j++) {
4613 fxOverX0ShieldTrks[j] = -1.;
4614 fxTimesRhoShieldTrks[j] = -1.;
4616 for(Int_t k=0; k<ntracks*6; k++) {
4617 fxOverX0LayerTrks[k] = -1.;
4618 fxTimesRhoLayerTrks[k] = -1.;
4625 //------------------------------------------------------------------------
4626 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4627 //-----------------------------------------------------------------
4628 // Delete LUT for storing material for each prolonged track
4629 //-----------------------------------------------------------------
4630 if(fxOverX0PipeTrks) {
4631 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4633 if(fxOverX0ShieldTrks) {
4634 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4637 if(fxOverX0LayerTrks) {
4638 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4640 if(fxTimesRhoPipeTrks) {
4641 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4643 if(fxTimesRhoShieldTrks) {
4644 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4646 if(fxTimesRhoLayerTrks) {
4647 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4651 //------------------------------------------------------------------------
4652 void AliITStrackerMI::SetForceSkippingOfLayer() {
4653 //-----------------------------------------------------------------
4654 // Check if we are forced to skip layers
4655 // either we set to skip them in RecoParam
4656 // or they were off during data-taking
4657 //-----------------------------------------------------------------
4659 const AliEventInfo *eventInfo = GetEventInfo();
4661 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4662 fForceSkippingOfLayer[l] = 0;
4664 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4668 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4669 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4671 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4672 } else if(l==2 || l==3) {
4673 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4675 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4681 //------------------------------------------------------------------------
4682 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4683 Int_t ilayer,Int_t idet) const {
4684 //-----------------------------------------------------------------
4685 // This method is used to decide whether to allow a prolongation
4686 // without clusters, because we want to skip the layer.
4687 // In this case the return value is > 0:
4688 // return 1: the user requested to skip a layer
4689 // return 2: track outside z acceptance
4690 //-----------------------------------------------------------------
4692 if (ForceSkippingOfLayer(ilayer)) return 1;
4694 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4696 if (idet<0 && // out in z
4697 ilayer>innerLayCanSkip &&
4698 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4699 // check if track will cross SPD outer layer
4700 Double_t phiAtSPD2,zAtSPD2;
4701 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4702 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4704 return 2; // always allow skipping, changed on 05.11.2009
4709 //------------------------------------------------------------------------
4710 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4711 Int_t ilayer,Int_t idet,
4712 Double_t dz,Double_t dy,
4713 Bool_t noClusters) const {
4714 //-----------------------------------------------------------------
4715 // This method is used to decide whether to allow a prolongation
4716 // without clusters, because there is a dead zone in the road.
4717 // In this case the return value is > 0:
4718 // return 1: dead zone at z=0,+-7cm in SPD
4719 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4720 // return 2: all road is "bad" (dead or noisy) from the OCDB
4721 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4722 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4723 //-----------------------------------------------------------------
4725 // check dead zones at z=0,+-7cm in the SPD
4726 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4727 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4728 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4729 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4730 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4731 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4732 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4733 for (Int_t i=0; i<3; i++)
4734 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4735 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4736 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4740 // check bad zones from OCDB
4741 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4743 if (idet<0) return 0;
4745 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4748 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4749 if (ilayer==0 || ilayer==1) { // ---------- SPD
4751 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4753 detSizeFactorX *= 2.;
4754 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4757 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4758 if (detType==2) segm->SetLayer(ilayer+1);
4759 Float_t detSizeX = detSizeFactorX*segm->Dx();
4760 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4762 // check if the road overlaps with bad chips
4764 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4765 Float_t zlocmin = zloc-dz;
4766 Float_t zlocmax = zloc+dz;
4767 Float_t xlocmin = xloc-dy;
4768 Float_t xlocmax = xloc+dy;
4769 Int_t chipsInRoad[100];
4771 // check if road goes out of detector
4772 Bool_t touchNeighbourDet=kFALSE;
4773 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4774 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4775 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4776 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4777 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));
4779 // check if this detector is bad
4781 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4782 if(!touchNeighbourDet) {
4783 return 2; // all detectors in road are bad
4785 return 3; // at least one is bad
4789 if(zlocmin>zlocmax)return 0;
4790 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4791 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4792 if (!nChipsInRoad) return 0;
4794 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4795 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4796 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4797 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4798 if (det.IsChipBad(chipsInRoad[iCh])) {
4806 if(!touchNeighbourDet) {
4807 AliDebug(2,"all bad in road");
4808 return 2; // all chips in road are bad
4810 return 3; // at least a bad chip in road
4815 AliDebug(2,"at least a bad in road");
4816 return 3; // at least a bad chip in road
4820 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4821 || !noClusters) return 0;
4823 // There are no clusters in road: check if there is at least
4824 // a bad SPD pixel or SDD anode or SSD strips on both sides
4826 Int_t idetInITS=idet;
4827 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4829 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4830 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4833 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4837 //------------------------------------------------------------------------
4838 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4839 const AliITStrackMI *track,
4840 Float_t &xloc,Float_t &zloc) const {
4841 //-----------------------------------------------------------------
4842 // Gives position of track in local module ref. frame
4843 //-----------------------------------------------------------------
4848 if(idet<0) return kTRUE; // track out of z acceptance of layer
4850 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4852 Int_t lad = Int_t(idet/ndet) + 1;
4854 Int_t det = idet - (lad-1)*ndet + 1;
4856 Double_t xyzGlob[3],xyzLoc[3];
4858 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4859 // take into account the misalignment: xyz at real detector plane
4860 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4862 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4864 xloc = (Float_t)xyzLoc[0];
4865 zloc = (Float_t)xyzLoc[2];
4869 //------------------------------------------------------------------------
4870 //------------------------------------------------------------------------
4871 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4873 // Method to be optimized further:
4874 // Aim: decide whether a track can be used for PlaneEff evaluation
4875 // the decision is taken based on the track quality at the layer under study
4876 // no information on the clusters on this layer has to be used
4877 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4878 // the cut is done on number of sigmas from the boundaries
4880 // Input: Actual track, layer [0,5] under study
4882 // Return: kTRUE if this is a good track
4884 // it will apply a pre-selection to obtain good quality tracks.
4885 // Here also you will have the possibility to put a control on the
4886 // impact point of the track on the basic block, in order to exclude border regions
4887 // this will be done by calling a proper method of the AliITSPlaneEff class.
4889 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4890 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4892 Int_t index[AliITSgeomTGeo::kNLayers];
4894 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4896 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4897 index[k]=clusters[k];
4901 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4902 AliITSlayer &layer=fgLayers[ilayer];
4903 Double_t r=layer.GetR();
4904 AliITStrackMI tmp(*track);
4906 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4907 Int_t nclout=0; Int_t nclin=0;
4908 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4909 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4910 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4911 // if (tmp.GetClIndex(lay)>=0) nclout++;
4912 if(index[lay]>=0)nclout++;
4914 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4915 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4916 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4917 if (index[lay]>=0) nclin++;
4919 Int_t ncl=nclout+nclin;
4920 Bool_t nextout = kFALSE;
4921 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4922 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4923 Bool_t nextin = kFALSE;
4924 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4925 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4926 // maximum number of missing clusters allowed in outermost layers
4927 if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4929 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4930 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4932 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4933 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4934 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4935 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4939 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4940 Int_t idet=layer.FindDetectorIndex(phi,z);
4941 if(idet<0) { AliInfo(Form("cannot find detector"));
4944 // here check if it has good Chi Square.
4946 //propagate to the intersection with the detector plane
4947 const AliITSdetector &det=layer.GetDetector(idet);
4948 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4952 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4953 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4954 if(key>fPlaneEff->Nblock()) return kFALSE;
4955 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4956 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4958 // DEFINITION OF SEARCH ROAD FOR accepting a track
4960 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4961 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4962 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4963 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4964 // done for RecPoints
4966 // exclude tracks at boundary between detectors
4967 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4968 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4969 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4970 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4971 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4972 if ( (locx-dx < blockXmn+boundaryWidth) ||
4973 (locx+dx > blockXmx-boundaryWidth) ||
4974 (locz-dz < blockZmn+boundaryWidth) ||
4975 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4978 //------------------------------------------------------------------------
4979 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4981 // This Method has to be optimized! For the time-being it uses the same criteria
4982 // as those used in the search of extra clusters for overlapping modules.
4984 // Method Purpose: estabilish whether a track has produced a recpoint or not
4985 // in the layer under study (For Plane efficiency)
4987 // inputs: AliITStrackMI* track (pointer to a usable track)
4989 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4990 // traversed by this very track. In details:
4991 // - if a cluster can be associated to the track then call
4992 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4993 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4996 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4997 AliITSlayer &layer=fgLayers[ilayer];
4998 Double_t r=layer.GetR();
4999 AliITStrackMI tmp(*track);
5003 if (!tmp.GetPhiZat(r,phi,z)) return;
5004 Int_t idet=layer.FindDetectorIndex(phi,z);
5006 if(idet<0) { AliInfo(Form("cannot find detector"));
5010 //propagate to the intersection with the detector plane
5011 const AliITSdetector &det=layer.GetDetector(idet);
5012 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5016 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5018 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5019 TMath::Sqrt(tmp.GetSigmaZ2() +
5020 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5021 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5022 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5023 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5024 TMath::Sqrt(tmp.GetSigmaY2() +
5025 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5026 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5027 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5029 // road in global (rphi,z) [i.e. in tracking ref. system]
5030 Double_t zmin = tmp.GetZ() - dz;
5031 Double_t zmax = tmp.GetZ() + dz;
5032 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5033 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5035 // select clusters in road
5036 layer.SelectClusters(zmin,zmax,ymin,ymax);
5038 // Define criteria for track-cluster association
5039 Double_t msz = tmp.GetSigmaZ2() +
5040 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5041 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5042 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5043 Double_t msy = tmp.GetSigmaY2() +
5044 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5045 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5046 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5047 if (tmp.GetConstrain()) {
5048 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5049 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5051 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5052 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5054 msz = 1./msz; // 1/RoadZ^2
5055 msy = 1./msy; // 1/RoadY^2
5058 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5060 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5061 //Double_t tolerance=0.2;
5062 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5063 idetc = cl->GetDetectorIndex();
5064 if(idet!=idetc) continue;
5065 //Int_t ilay = cl->GetLayer();
5067 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5068 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5070 Double_t chi2=tmp.GetPredictedChi2(cl);
5071 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5075 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5077 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5078 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5079 if(key>fPlaneEff->Nblock()) return;
5080 Bool_t found=kFALSE;
5083 while ((cl=layer.GetNextCluster(clidx))!=0) {
5084 idetc = cl->GetDetectorIndex();
5085 if(idet!=idetc) continue;
5086 // here real control to see whether the cluster can be associated to the track.
5087 // cluster not associated to track
5088 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5089 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5090 // calculate track-clusters chi2
5091 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5092 // in particular, the error associated to the cluster
5093 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5095 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5097 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5098 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5099 // track->SetExtraModule(ilayer,idetExtra);
5101 if(!fPlaneEff->UpDatePlaneEff(found,key))
5102 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5104 // this for FO efficiency studies (only for SPD) //
5105 UInt_t keyFO=999999;
5106 Bool_t foundFO=kFALSE;
5107 if(ilayer<2){ //ONLY SPD layers for FastOr studies
5108 TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5109 Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5110 if(!fSPDChipIntPlaneEff[key]){
5111 AliITSPlaneEffSPD spd;
5112 keyFO = spd.SwitchChipKeyNumbering(key);
5113 if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5114 keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5115 if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5116 AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5119 if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5120 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5126 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5127 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5128 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5129 Int_t cltype[2]={-999,-999};
5132 Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5136 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5137 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5140 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5141 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5142 cltype[0]=layer.GetCluster(ci)->GetNy();
5143 cltype[1]=layer.GetCluster(ci)->GetNz();
5145 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5146 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5147 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5148 // It is computed properly by calling the method
5149 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5151 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5152 //if (tmp.PropagateTo(x,0.,0.)) {
5153 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5154 AliCluster c(*layer.GetCluster(ci));
5155 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5156 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5157 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5158 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5159 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5162 // Compute the angles between the track and the module
5163 // compute the angle "in phi direction", i.e. the angle in the transverse plane
5164 // between the normal to the module and the projection (in the transverse plane) of the
5166 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5167 Float_t tgl = tmp.GetTgl();
5168 Float_t phitr = tmp.GetSnp();
5169 phitr = TMath::ASin(phitr);
5170 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5172 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5173 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5175 alpha = tmp.GetAlpha();
5176 Double_t phiglob = alpha+phitr;
5178 p[0] = TMath::Cos(phiglob);
5179 p[1] = TMath::Sin(phiglob);
5181 TVector3 pvec(p[0],p[1],p[2]);
5182 TVector3 normvec(rot[1],rot[4],rot[7]);
5183 Double_t angle = pvec.Angle(normvec);
5185 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5186 angle *= 180./TMath::Pi();
5189 TVector3 pt(p[0],p[1],0);
5190 TVector3 normt(rot[1],rot[4],0);
5191 Double_t anglet = pt.Angle(normt);
5193 Double_t phiPt = TMath::ATan2(p[1],p[0]);
5194 if(phiPt<0)phiPt+=2.*TMath::Pi();
5195 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5196 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5197 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5198 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
5199 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5200 anglet *= 180./TMath::Pi();
5202 angleModTrack[2]=(Float_t) angle;
5203 angleModTrack[0]=(Float_t) anglet;
5204 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5205 angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5206 angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5207 angleModTrack[1]*=180./TMath::Pi(); // in degree
5209 fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
5211 // For FO efficiency studies of SPD
5212 if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
5214 if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5218 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5220 // find the MC cluster for the label
5221 return fgLayers[lr].FindClusterForLabel(label,store);
5225 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5227 // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5228 strncpy(patt,"......",6);
5230 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5233 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5234 Int_t cindex = track->GetClusterIndex(i);
5235 Int_t l=(cindex & 0xf0000000) >> 28;
5236 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5238 for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5239 patt[l] = isWrong ? 'f':'c';
5245 //------------------------------------------------------------------------
5246 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5248 //--------------------------------------------------------------------
5251 for (int ic=0;ic<fN;ic++) {
5252 const AliITSRecPoint *cl = GetCluster(ic);
5253 for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5255 if (store) store[nfound] = ic;