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);
710 double bz0 = GetBz();
711 const double kWatchStep=10.; // for larger steps watch arc vs segment difference
714 for (Int_t i=0; i<nentr; i++) {
715 AliESDtrack *esd=event->GetTrack(i);
717 // Start time integral and add distance from current position to vertex
718 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
719 AliITStrackMI t(*esd);
720 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
724 double dxs = xyzTrk[0] - xyzVtx[0];
725 double dys = xyzTrk[1] - xyzVtx[1];
726 double dzs = xyzTrk[2] - xyzVtx[2];
727 // RS: for large segment steps d use approximation of cicrular arc s by
728 // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
729 // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
730 // Hence s^2/d^2 = (1+1/6 p^2)^2
731 dst2 = dxs*dxs + dys*dys;
732 if (dst2 > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
733 double crv = TMath::Abs(esd->GetC(bz0));
734 double fcarc = 1.+crv*crv*dst2/6.;
739 t.StartTimeIntegral();
740 t.AddTimeStep(TMath::Sqrt(dst2));
742 // transfer the time integral to ESD track
743 esd->SetStatus(AliESDtrack::kTIME);
744 Double_t times[10];t.GetIntegratedTimes(times); esd->SetIntegratedTimes(times);
745 esd->SetIntegratedLength(t.GetIntegratedLength());
747 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
749 t.SetExpQ(TMath::Max(0.8*t.GetESDtrack()->GetTPCsignal(),30.));
750 ResetTrackToFollow(t);
752 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
753 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,&t)) {
754 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) continue;
755 fTrackToFollow.SetLabel(t.GetLabel());
756 //fTrackToFollow.CookdEdx();
757 CookLabel(&fTrackToFollow,0.); //For comparison only
758 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
759 //UseClusters(&fTrackToFollow);
764 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
766 fTrackingPhase="Default";
770 //------------------------------------------------------------------------
771 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
772 //--------------------------------------------------------------------
773 // This functions refits ITS tracks using the
774 // "inward propagated" TPC tracks
775 // The clusters must be loaded !
776 //--------------------------------------------------------------------
777 fTrackingPhase="RefitInward";
779 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
781 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
782 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
784 Int_t nentr=event->GetNumberOfTracks();
785 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
787 // only for PlaneEff and in case of SPD (for FO studies)
788 if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
789 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 &&
790 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) {
791 for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
795 for (Int_t i=0; i<nentr; i++) {
796 AliESDtrack *esd=event->GetTrack(i);
798 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
799 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
800 if (esd->GetStatus()&AliESDtrack::kTPCout)
801 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
803 AliITStrackMI *t = new AliITStrackMI(*esd);
805 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
806 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
811 ResetTrackToFollow(*t);
812 fTrackToFollow.ResetClusters();
814 // ITS standalone tracks
815 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
816 fTrackToFollow.ResetCovariance(10.);
817 // protection for loopers that can have parameters screwed up
818 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
819 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
826 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
827 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
829 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
830 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
831 AliDebug(2," refit OK");
832 fTrackToFollow.SetLabel(t->GetLabel());
833 // fTrackToFollow.CookdEdx();
834 CookdEdx(&fTrackToFollow);
836 CookLabel(&fTrackToFollow,0.0); //For comparison only
839 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
840 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
841 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
842 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
843 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
844 Double_t r[3]={0.,0.,0.};
846 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
853 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
855 fTrackingPhase="Default";
859 //------------------------------------------------------------------------
860 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
861 //--------------------------------------------------------------------
862 // Return pointer to a given cluster
863 //--------------------------------------------------------------------
864 Int_t l=(index & 0xf0000000) >> 28;
865 Int_t c=(index & 0x0fffffff) >> 00;
866 return fgLayers[l].GetCluster(c);
868 //------------------------------------------------------------------------
869 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
870 //--------------------------------------------------------------------
871 // Get track space point with index i
872 //--------------------------------------------------------------------
874 Int_t l=(index & 0xf0000000) >> 28;
875 Int_t c=(index & 0x0fffffff) >> 00;
876 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
877 Int_t idet = cl->GetDetectorIndex();
881 cl->GetGlobalXYZ(xyz);
882 cl->GetGlobalCov(cov);
884 p.SetCharge(cl->GetQ());
885 p.SetDriftTime(cl->GetDriftTime());
886 p.SetChargeRatio(cl->GetChargeRatio());
887 p.SetClusterType(cl->GetClusterType());
888 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
891 iLayer = AliGeomManager::kSPD1;
894 iLayer = AliGeomManager::kSPD2;
897 iLayer = AliGeomManager::kSDD1;
900 iLayer = AliGeomManager::kSDD2;
903 iLayer = AliGeomManager::kSSD1;
906 iLayer = AliGeomManager::kSSD2;
909 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
912 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
913 p.SetVolumeID((UShort_t)volid);
916 //------------------------------------------------------------------------
917 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
918 AliTrackPoint& p, const AliESDtrack *t) {
919 //--------------------------------------------------------------------
920 // Get track space point with index i
921 // (assign error estimated during the tracking)
922 //--------------------------------------------------------------------
924 Int_t l=(index & 0xf0000000) >> 28;
925 Int_t c=(index & 0x0fffffff) >> 00;
926 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
927 Int_t idet = cl->GetDetectorIndex();
929 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
931 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
933 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
934 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
935 Double_t alpha = t->GetAlpha();
936 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
937 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
938 phi += alpha-det.GetPhi();
939 Float_t tgphi = TMath::Tan(phi);
941 Float_t tgl = t->GetTgl(); // tgl about const along track
942 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
944 Float_t errtrky,errtrkz,covyz;
945 Bool_t addMisalErr=kFALSE;
946 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
950 cl->GetGlobalXYZ(xyz);
951 // cl->GetGlobalCov(cov);
952 Float_t pos[3] = {0.,0.,0.};
953 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
954 tmpcl.GetGlobalCov(cov);
957 p.SetCharge(cl->GetQ());
958 p.SetDriftTime(cl->GetDriftTime());
959 p.SetChargeRatio(cl->GetChargeRatio());
960 p.SetClusterType(cl->GetClusterType());
962 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
965 iLayer = AliGeomManager::kSPD1;
968 iLayer = AliGeomManager::kSPD2;
971 iLayer = AliGeomManager::kSDD1;
974 iLayer = AliGeomManager::kSDD2;
977 iLayer = AliGeomManager::kSSD1;
980 iLayer = AliGeomManager::kSSD2;
983 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
986 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
988 p.SetVolumeID((UShort_t)volid);
991 //------------------------------------------------------------------------
992 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
994 //--------------------------------------------------------------------
995 // Follow prolongation tree
996 //--------------------------------------------------------------------
998 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
999 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
1002 AliESDtrack * esd = otrack->GetESDtrack();
1003 if (esd->GetV0Index(0)>0) {
1004 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
1005 // mapping of ESD track is different as ITS track in Containers
1006 // Need something more stable
1007 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
1008 for (Int_t i=0;i<3;i++){
1009 Int_t index = esd->GetV0Index(i);
1010 if (index==0) break;
1011 AliESDv0 * vertex = fEsd->GetV0(index);
1012 if (vertex->GetStatus()<0) continue; // rejected V0
1014 if (esd->GetSign()>0) {
1015 vertex->SetIndex(0,esdindex);
1017 vertex->SetIndex(1,esdindex);
1021 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1023 bestarray = new TObjArray(5);
1024 bestarray->SetOwner();
1025 fBestHypothesys.AddAt(bestarray,esdindex);
1029 //setup tree of the prolongations
1031 const int kMaxTr = 100; //RS
1032 static AliITStrackMI tracks[7][kMaxTr];
1033 AliITStrackMI *currenttrack;
1034 static AliITStrackMI currenttrack1;
1035 static AliITStrackMI currenttrack2;
1036 static AliITStrackMI backuptrack;
1038 Int_t nindexes[7][kMaxTr];
1039 Float_t normalizedchi2[kMaxTr];
1040 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1041 otrack->SetNSkipped(0);
1042 new (&(tracks[6][0])) AliITStrackMI(*otrack);
1044 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1045 Int_t modstatus = 1; // found
1049 // follow prolongations
1050 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1051 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1054 AliITSlayer &layer=fgLayers[ilayer];
1055 Double_t r = layer.GetR();
1061 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1062 // printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1064 if (ntracks[ilayer]>=kMaxTr) break;
1065 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1066 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1067 if (ntracks[ilayer]>15+ilayer){
1068 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1069 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1072 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1074 // material between SSD and SDD, SDD and SPD
1076 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1078 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1082 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1083 Int_t idet=layer.FindDetectorIndex(phi,z);
1085 Double_t trackGlobXYZ1[3];
1086 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1088 // Get the budget to the primary vertex for the current track being prolonged
1089 Double_t budgetToPrimVertex = 0;
1090 double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1093 if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1094 else budgetToPrimVertex = GetEffectiveThickness();
1096 // check if we allow a prolongation without point
1097 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1099 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1100 // propagate to the layer radius
1101 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1102 if(!vtrack->Propagate(xToGo)) continue;
1103 // apply correction for material of the current layer
1104 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1105 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1106 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1107 vtrack->SetClIndex(ilayer,-1);
1108 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1109 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1110 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1112 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1114 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1115 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1121 // track outside layer acceptance in z
1122 if (idet<0) continue;
1124 //propagate to the intersection with the detector plane
1125 const AliITSdetector &det=layer.GetDetector(idet);
1126 new(¤ttrack2) AliITStrackMI(currenttrack1);
1127 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1128 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1129 currenttrack1.SetDetectorIndex(idet);
1130 currenttrack2.SetDetectorIndex(idet);
1131 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1134 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1136 // road in global (rphi,z) [i.e. in tracking ref. system]
1137 Double_t zmin,zmax,ymin,ymax;
1138 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1140 // select clusters in road
1141 layer.SelectClusters(zmin,zmax,ymin,ymax);
1142 //********************
1144 // Define criteria for track-cluster association
1145 Double_t msz = currenttrack1.GetSigmaZ2() +
1146 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1147 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1148 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1149 Double_t msy = currenttrack1.GetSigmaY2() +
1150 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1151 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1152 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1154 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1155 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1157 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1158 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1160 msz = 1./msz; // 1/RoadZ^2
1161 msy = 1./msy; // 1/RoadY^2
1165 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1167 const AliITSRecPoint *cl=0;
1169 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1170 Bool_t deadzoneSPD=kFALSE;
1171 currenttrack = ¤ttrack1;
1173 // check if the road contains a dead zone
1174 Bool_t noClusters = kFALSE;
1175 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1176 if (noClusters) AliDebug(2,"no clusters in road");
1177 Double_t dz=0.5*(zmax-zmin);
1178 Double_t dy=0.5*(ymax-ymin);
1179 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1180 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1181 // create a prolongation without clusters (check also if there are no clusters in the road)
1184 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1185 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1186 updatetrack->SetClIndex(ilayer,-1);
1188 modstatus = 5; // no cls in road
1189 } else if (dead==1) {
1190 modstatus = 7; // holes in z in SPD
1191 } else if (dead==2 || dead==3 || dead==4) {
1192 modstatus = 2; // dead from OCDB
1194 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1195 // apply correction for material of the current layer
1196 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1197 if (constrain) { // apply vertex constrain
1198 updatetrack->SetConstrain(constrain);
1199 Bool_t isPrim = kTRUE;
1200 if (ilayer<4) { // check that it's close to the vertex
1201 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1202 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1203 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1204 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1205 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1207 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1209 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1210 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1213 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1215 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1216 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1218 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1219 updatetrack->SetDeadZoneProbability(ilayer,1.);
1220 } else if (dead==4) { // at least a single dead channel from OCDB
1221 updatetrack->SetDeadZoneProbability(ilayer,0.);
1228 // loop over clusters in the road
1229 while ((cl=layer.GetNextCluster(clidx))!=0) {
1230 if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters
1231 Bool_t changedet =kFALSE;
1232 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1233 Int_t idetc=cl->GetDetectorIndex();
1235 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1236 // take into account misalignment (bring track to real detector plane)
1237 Double_t xTrOrig = currenttrack->GetX();
1238 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1239 // a first cut on track-cluster distance
1240 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1241 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1242 { // cluster not associated to track
1243 AliDebug(2,"not associated");
1244 // MvL: added here as well
1245 // bring track back to ideal detector plane
1246 currenttrack->Propagate(xTrOrig);
1249 // bring track back to ideal detector plane
1250 if (!currenttrack->Propagate(xTrOrig)) continue;
1251 } else { // have to move track to cluster's detector
1252 const AliITSdetector &detc=layer.GetDetector(idetc);
1253 // a first cut on track-cluster distance
1255 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1256 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1257 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1258 continue; // cluster not associated to track
1260 new (&backuptrack) AliITStrackMI(currenttrack2);
1262 currenttrack =¤ttrack2;
1263 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1264 new (currenttrack) AliITStrackMI(backuptrack);
1268 currenttrack->SetDetectorIndex(idetc);
1269 // Get again the budget to the primary vertex
1270 // for the current track being prolonged, if had to change detector
1271 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1274 // calculate track-clusters chi2
1275 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1277 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1278 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1279 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1280 if (ntracks[ilayer]>=kMaxTr) continue;
1281 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1282 updatetrack->SetClIndex(ilayer,-1);
1283 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1285 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1286 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1287 AliDebug(2,"update failed");
1290 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1291 modstatus = 1; // found
1292 } else { // virtual cluster in dead zone
1293 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1294 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1295 modstatus = 7; // holes in z in SPD
1299 Float_t xlocnewdet,zlocnewdet;
1300 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1301 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1304 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1306 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1308 // apply correction for material of the current layer
1309 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1311 if (constrain) { // apply vertex constrain
1312 updatetrack->SetConstrain(constrain);
1313 Bool_t isPrim = kTRUE;
1314 if (ilayer<4) { // check that it's close to the vertex
1315 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1316 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1317 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1318 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1319 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1321 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1323 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1324 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1326 } //apply vertex constrain
1328 } // create new hypothesis
1330 AliDebug(2,"chi2 too large");
1333 } // loop over possible prolongations
1335 // allow one prolongation without clusters
1336 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1337 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1338 // apply correction for material of the current layer
1339 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1340 vtrack->SetClIndex(ilayer,-1);
1341 modstatus = 3; // skipped
1342 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1343 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1345 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1346 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1348 vtrack->IncrementNSkipped();
1353 } // loop over tracks in layer ilayer+1
1355 //loop over track candidates for the current layer
1361 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1362 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1363 if (normalizedchi2[itrack] <
1364 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1368 if (constrain) { // constrain
1369 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1371 } else { // !constrain
1372 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1377 // sort tracks by increasing normalized chi2
1378 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1379 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1380 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1381 // if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1382 if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9);
1383 } // end loop over layers
1387 // Now select tracks to be kept
1389 Int_t max = constrain ? 20 : 5;
1391 // tracks that reach layer 0 (SPD inner)
1392 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1393 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1394 if (track.GetNumberOfClusters()<2) continue;
1395 if (!constrain && track.GetNormChi2(0) >
1396 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1399 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1402 // tracks that reach layer 1 (SPD outer)
1403 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1404 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1405 if (track.GetNumberOfClusters()<4) continue;
1406 if (!constrain && track.GetNormChi2(1) >
1407 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1408 if (constrain) track.IncrementNSkipped();
1410 track.SetD(0,track.GetD(GetX(),GetY()));
1411 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1412 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1413 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1416 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1419 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1421 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1422 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1423 if (track.GetNumberOfClusters()<3) continue;
1424 if (track.GetNormChi2(2) >
1425 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1426 track.SetD(0,track.GetD(GetX(),GetY()));
1427 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1428 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1429 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1431 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1437 // register best track of each layer - important for V0 finder
1439 for (Int_t ilayer=0;ilayer<5;ilayer++){
1440 if (ntracks[ilayer]==0) continue;
1441 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1442 if (track.GetNumberOfClusters()<1) continue;
1443 CookLabel(&track,0);
1444 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1448 // update TPC V0 information
1450 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1451 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1452 for (Int_t i=0;i<3;i++){
1453 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1454 if (index==0) break;
1455 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1456 if (vertex->GetStatus()<0) continue; // rejected V0
1458 if (otrack->GetSign()>0) {
1459 vertex->SetIndex(0,esdindex);
1462 vertex->SetIndex(1,esdindex);
1464 //find nearest layer with track info
1465 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1466 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1467 Int_t nearest = nearestold;
1468 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1469 if (ntracks[nearest]==0){
1474 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1475 if (nearestold<5&&nearest<5){
1476 Bool_t accept = track.GetNormChi2(nearest)<10;
1478 if (track.GetSign()>0) {
1479 vertex->SetParamP(track);
1480 vertex->Update(fprimvertex);
1481 //vertex->SetIndex(0,track.fESDtrack->GetID());
1482 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1484 vertex->SetParamN(track);
1485 vertex->Update(fprimvertex);
1486 //vertex->SetIndex(1,track.fESDtrack->GetID());
1487 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1489 vertex->SetStatus(vertex->GetStatus()+1);
1491 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1498 //------------------------------------------------------------------------
1499 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1501 //--------------------------------------------------------------------
1504 return fgLayers[layer];
1506 //------------------------------------------------------------------------
1507 AliITStrackerMI::AliITSlayer::AliITSlayer():
1537 //--------------------------------------------------------------------
1538 //default AliITSlayer constructor
1539 //--------------------------------------------------------------------
1540 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1541 fClusterWeight[i]=0;
1542 fClusterTracks[0][i]=-1;
1543 fClusterTracks[1][i]=-1;
1544 fClusterTracks[2][i]=-1;
1545 fClusterTracks[3][i]=-1;
1552 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1553 for (Int_t j1=0; j1<6; j1++) {
1554 fClusters5[j1][j]=0;
1555 fClusterIndex5[j1][j]=-1;
1564 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1565 for (Int_t j1=0; j1<11; j1++) {
1566 fClusters10[j1][j]=0;
1567 fClusterIndex10[j1][j]=-1;
1576 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1577 for (Int_t j1=0; j1<21; j1++) {
1578 fClusters20[j1][j]=0;
1579 fClusterIndex20[j1][j]=-1;
1587 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1592 //------------------------------------------------------------------------
1593 AliITStrackerMI::AliITSlayer::
1594 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1623 //--------------------------------------------------------------------
1624 //main AliITSlayer constructor
1625 //--------------------------------------------------------------------
1626 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1627 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1629 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1630 fClusterWeight[i]=0;
1631 fClusterTracks[0][i]=-1;
1632 fClusterTracks[1][i]=-1;
1633 fClusterTracks[2][i]=-1;
1634 fClusterTracks[3][i]=-1;
1642 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1643 for (Int_t j1=0; j1<6; j1++) {
1644 fClusters5[j1][j]=0;
1645 fClusterIndex5[j1][j]=-1;
1654 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1655 for (Int_t j1=0; j1<11; j1++) {
1656 fClusters10[j1][j]=0;
1657 fClusterIndex10[j1][j]=-1;
1666 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1667 for (Int_t j1=0; j1<21; j1++) {
1668 fClusters20[j1][j]=0;
1669 fClusterIndex20[j1][j]=-1;
1677 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1683 //------------------------------------------------------------------------
1684 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1686 fPhiOffset(layer.fPhiOffset),
1687 fNladders(layer.fNladders),
1688 fZOffset(layer.fZOffset),
1689 fNdetectors(layer.fNdetectors),
1690 fDetectors(layer.fDetectors),
1695 fClustersCs(layer.fClustersCs),
1696 fClusterIndexCs(layer.fClusterIndexCs),
1700 fCurrentSlice(layer.fCurrentSlice),
1708 fAccepted(layer.fAccepted),
1710 fMaxSigmaClY(layer.fMaxSigmaClY),
1711 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1712 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1717 //------------------------------------------------------------------------
1718 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1719 //--------------------------------------------------------------------
1720 // AliITSlayer destructor
1721 //--------------------------------------------------------------------
1722 delete [] fDetectors;
1723 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1724 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1725 fClusterWeight[i]=0;
1726 fClusterTracks[0][i]=-1;
1727 fClusterTracks[1][i]=-1;
1728 fClusterTracks[2][i]=-1;
1729 fClusterTracks[3][i]=-1;
1732 //------------------------------------------------------------------------
1733 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1734 //--------------------------------------------------------------------
1735 // This function removes loaded clusters
1736 //--------------------------------------------------------------------
1737 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1738 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1739 fClusterWeight[i]=0;
1740 fClusterTracks[0][i]=-1;
1741 fClusterTracks[1][i]=-1;
1742 fClusterTracks[2][i]=-1;
1743 fClusterTracks[3][i]=-1;
1749 //------------------------------------------------------------------------
1750 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1751 //--------------------------------------------------------------------
1752 // This function reset weights of the clusters
1753 //--------------------------------------------------------------------
1754 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1755 fClusterWeight[i]=0;
1756 fClusterTracks[0][i]=-1;
1757 fClusterTracks[1][i]=-1;
1758 fClusterTracks[2][i]=-1;
1759 fClusterTracks[3][i]=-1;
1761 for (Int_t i=0; i<fN;i++) {
1762 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1763 if (cl&&cl->IsUsed()) cl->Use();
1767 //------------------------------------------------------------------------
1768 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1769 //--------------------------------------------------------------------
1770 // This function calculates the road defined by the cluster density
1771 //--------------------------------------------------------------------
1773 for (Int_t i=0; i<fN; i++) {
1774 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1776 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1778 //------------------------------------------------------------------------
1779 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1780 //--------------------------------------------------------------------
1781 //This function adds a cluster to this layer
1782 //--------------------------------------------------------------------
1783 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1789 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1791 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1792 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1793 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1794 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1795 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1796 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1799 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1800 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1801 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1802 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1806 //------------------------------------------------------------------------
1807 void AliITStrackerMI::AliITSlayer::SortClusters()
1812 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1813 Float_t *z = new Float_t[fN];
1814 Int_t * index = new Int_t[fN];
1816 fMaxSigmaClY=0.; //AD
1817 fMaxSigmaClZ=0.; //AD
1819 for (Int_t i=0;i<fN;i++){
1820 z[i] = fClusters[i]->GetZ();
1821 // save largest errors in y and z for this layer
1822 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1823 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1825 TMath::Sort(fN,z,index,kFALSE);
1826 for (Int_t i=0;i<fN;i++){
1827 clusters[i] = fClusters[index[i]];
1830 for (Int_t i=0;i<fN;i++){
1831 fClusters[i] = clusters[i];
1832 fZ[i] = fClusters[i]->GetZ();
1833 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1834 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1835 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1845 for (Int_t i=0;i<fN;i++){
1846 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1847 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1848 fClusterIndex[i] = i;
1852 fDy5 = (fYB[1]-fYB[0])/5.;
1853 fDy10 = (fYB[1]-fYB[0])/10.;
1854 fDy20 = (fYB[1]-fYB[0])/20.;
1855 for (Int_t i=0;i<6;i++) fN5[i] =0;
1856 for (Int_t i=0;i<11;i++) fN10[i]=0;
1857 for (Int_t i=0;i<21;i++) fN20[i]=0;
1859 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;}
1860 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;}
1861 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;}
1864 for (Int_t i=0;i<fN;i++)
1865 for (Int_t irot=-1;irot<=1;irot++){
1866 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1868 for (Int_t slice=0; slice<6;slice++){
1869 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1870 fClusters5[slice][fN5[slice]] = fClusters[i];
1871 fY5[slice][fN5[slice]] = curY;
1872 fZ5[slice][fN5[slice]] = fZ[i];
1873 fClusterIndex5[slice][fN5[slice]]=i;
1878 for (Int_t slice=0; slice<11;slice++){
1879 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1880 fClusters10[slice][fN10[slice]] = fClusters[i];
1881 fY10[slice][fN10[slice]] = curY;
1882 fZ10[slice][fN10[slice]] = fZ[i];
1883 fClusterIndex10[slice][fN10[slice]]=i;
1888 for (Int_t slice=0; slice<21;slice++){
1889 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1890 fClusters20[slice][fN20[slice]] = fClusters[i];
1891 fY20[slice][fN20[slice]] = curY;
1892 fZ20[slice][fN20[slice]] = fZ[i];
1893 fClusterIndex20[slice][fN20[slice]]=i;
1900 // consistency check
1902 for (Int_t i=0;i<fN-1;i++){
1908 for (Int_t slice=0;slice<21;slice++)
1909 for (Int_t i=0;i<fN20[slice]-1;i++){
1910 if (fZ20[slice][i]>fZ20[slice][i+1]){
1917 //------------------------------------------------------------------------
1918 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1919 //--------------------------------------------------------------------
1920 // This function returns the index of the nearest cluster
1921 //--------------------------------------------------------------------
1924 if (fCurrentSlice<0) {
1933 if (ncl==0) return 0;
1934 Int_t b=0, e=ncl-1, m=(b+e)/2;
1935 for (; b<e; m=(b+e)/2) {
1936 // if (z > fClusters[m]->GetZ()) b=m+1;
1937 if (z > zcl[m]) b=m+1;
1942 //------------------------------------------------------------------------
1943 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 {
1944 //--------------------------------------------------------------------
1945 // This function computes the rectangular road for this track
1946 //--------------------------------------------------------------------
1949 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1950 // take into account the misalignment: propagate track to misaligned detector plane
1951 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1953 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1954 TMath::Sqrt(track->GetSigmaZ2() +
1955 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1956 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1957 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1958 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1959 TMath::Sqrt(track->GetSigmaY2() +
1960 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1961 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1962 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1964 // track at boundary between detectors, enlarge road
1965 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1966 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1967 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1968 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1969 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1970 Float_t tgl = TMath::Abs(track->GetTgl());
1971 if (tgl > 1.) tgl=1.;
1972 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1973 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1974 Float_t snp = TMath::Abs(track->GetSnp());
1975 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1976 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1979 // add to the road a term (up to 2-3 mm) to deal with misalignments
1980 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1981 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1983 Double_t r = fgLayers[ilayer].GetR();
1984 zmin = track->GetZ() - dz;
1985 zmax = track->GetZ() + dz;
1986 ymin = track->GetY() + r*det.GetPhi() - dy;
1987 ymax = track->GetY() + r*det.GetPhi() + dy;
1989 // bring track back to idead detector plane
1990 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1994 //------------------------------------------------------------------------
1995 void AliITStrackerMI::AliITSlayer::
1996 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1997 //--------------------------------------------------------------------
1998 // This function sets the "window"
1999 //--------------------------------------------------------------------
2001 Double_t circle=2*TMath::Pi()*fR;
2007 // enlarge road in y by maximum cluster error on this layer (3 sigma)
2008 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
2009 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
2010 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
2011 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
2013 Float_t ymiddle = (fYmin+fYmax)*0.5;
2014 if (ymiddle<fYB[0]) {
2015 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
2016 } else if (ymiddle>fYB[1]) {
2017 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
2023 fClustersCs = fClusters;
2024 fClusterIndexCs = fClusterIndex;
2030 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2031 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2032 if (slice<0) slice=0;
2033 if (slice>20) slice=20;
2034 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2036 fCurrentSlice=slice;
2037 fClustersCs = fClusters20[fCurrentSlice];
2038 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2039 fYcs = fY20[fCurrentSlice];
2040 fZcs = fZ20[fCurrentSlice];
2041 fNcs = fN20[fCurrentSlice];
2046 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2047 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2048 if (slice<0) slice=0;
2049 if (slice>10) slice=10;
2050 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2052 fCurrentSlice=slice;
2053 fClustersCs = fClusters10[fCurrentSlice];
2054 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2055 fYcs = fY10[fCurrentSlice];
2056 fZcs = fZ10[fCurrentSlice];
2057 fNcs = fN10[fCurrentSlice];
2062 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2063 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2064 if (slice<0) slice=0;
2065 if (slice>5) slice=5;
2066 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2068 fCurrentSlice=slice;
2069 fClustersCs = fClusters5[fCurrentSlice];
2070 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2071 fYcs = fY5[fCurrentSlice];
2072 fZcs = fZ5[fCurrentSlice];
2073 fNcs = fN5[fCurrentSlice];
2077 fI = FindClusterIndex(fZmin);
2078 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2084 //------------------------------------------------------------------------
2085 Int_t AliITStrackerMI::AliITSlayer::
2086 FindDetectorIndex(Double_t phi, Double_t z) const {
2087 //--------------------------------------------------------------------
2088 //This function finds the detector crossed by the track
2089 //--------------------------------------------------------------------
2091 if (fZOffset<0) // old geometry
2092 dphi = -(phi-fPhiOffset);
2093 else // new geometry
2094 dphi = phi-fPhiOffset;
2097 if (dphi < 0) dphi += 2*TMath::Pi();
2098 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2099 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2100 if (np>=fNladders) np-=fNladders;
2101 if (np<0) np+=fNladders;
2104 Double_t dz=fZOffset-z;
2105 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2106 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2107 if (nz>=fNdetectors || nz<0) {
2108 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2112 // ad hoc correction for 3rd ladder of SDD inner layer,
2113 // which is reversed (rotated by pi around local y)
2114 // this correction is OK only from AliITSv11Hybrid onwards
2115 if (GetR()>12. && GetR()<20.) { // SDD inner
2116 if(np==2) { // 3rd ladder
2117 Double_t posMod252[3];
2118 AliITSgeomTGeo::GetTranslation(252,posMod252);
2119 // check the Z coordinate of Mod 252: if negative
2120 // (old SDD geometry in AliITSv11Hybrid)
2121 // the swap of numeration whould be applied
2122 if(posMod252[2]<0.){
2123 nz = (fNdetectors-1) - nz;
2127 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2130 return np*fNdetectors + nz;
2132 //------------------------------------------------------------------------
2133 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2135 //--------------------------------------------------------------------
2136 // This function returns clusters within the "window"
2137 //--------------------------------------------------------------------
2139 if (fCurrentSlice<0) {
2140 Double_t rpi2 = 2.*fR*TMath::Pi();
2141 for (Int_t i=fI; i<fImax; i++) {
2144 if (fYmax<y) y -= rpi2;
2145 if (fYmin>y) y += rpi2;
2146 if (y<fYmin) continue;
2147 if (y>fYmax) continue;
2149 // skip clusters that are in "extended" road but they
2150 // 3sigma error does not touch the original road
2151 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2152 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2154 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2157 return fClusters[i];
2160 for (Int_t i=fI; i<fImax; i++) {
2161 if (fYcs[i]<fYmin) continue;
2162 if (fYcs[i]>fYmax) continue;
2163 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2164 ci=fClusterIndexCs[i];
2166 return fClustersCs[i];
2171 //------------------------------------------------------------------------
2172 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2174 //--------------------------------------------------------------------
2175 // This function returns the layer thickness at this point (units X0)
2176 //--------------------------------------------------------------------
2178 x0=AliITSRecoParam::GetX0Air();
2179 if (43<fR&&fR<45) { //SSD2
2182 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2183 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2184 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2185 for (Int_t i=0; i<12; i++) {
2186 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2187 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2191 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2192 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2196 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2197 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2200 if (37<fR&&fR<41) { //SSD1
2203 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2204 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2205 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2206 for (Int_t i=0; i<11; i++) {
2207 if (TMath::Abs(z-3.9*i)<0.15) {
2208 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2212 if (TMath::Abs(z+3.9*i)<0.15) {
2213 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2217 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2218 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2221 if (13<fR&&fR<26) { //SDD
2224 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2226 if (TMath::Abs(y-1.80)<0.55) {
2228 for (Int_t j=0; j<20; j++) {
2229 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2230 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2233 if (TMath::Abs(y+1.80)<0.55) {
2235 for (Int_t j=0; j<20; j++) {
2236 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2237 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2241 for (Int_t i=0; i<4; i++) {
2242 if (TMath::Abs(z-7.3*i)<0.60) {
2244 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2247 if (TMath::Abs(z+7.3*i)<0.60) {
2249 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2254 if (6<fR&&fR<8) { //SPD2
2255 Double_t dd=0.0063; x0=21.5;
2257 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2258 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2260 if (3<fR&&fR<5) { //SPD1
2261 Double_t dd=0.0063; x0=21.5;
2263 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2264 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2269 //------------------------------------------------------------------------
2270 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2272 fRmisal(det.fRmisal),
2274 fSinPhi(det.fSinPhi),
2275 fCosPhi(det.fCosPhi),
2281 fNChips(det.fNChips),
2282 fChipIsBad(det.fChipIsBad)
2286 //------------------------------------------------------------------------
2287 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2288 const AliITSDetTypeRec *detTypeRec)
2290 //--------------------------------------------------------------------
2291 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2292 //--------------------------------------------------------------------
2294 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2295 // while in the tracker they start from 0 for each layer
2296 for(Int_t il=0; il<ilayer; il++)
2297 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2300 if (ilayer==0 || ilayer==1) { // ---------- SPD
2302 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2304 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2307 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2311 // Get calibration from AliITSDetTypeRec
2312 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2313 calib->SetModuleIndex(idet);
2314 AliITSCalibration *calibSPDdead = 0;
2315 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2316 if (calib->IsBad() ||
2317 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2320 // printf("lay %d bad %d\n",ilayer,idet);
2323 // Get segmentation from AliITSDetTypeRec
2324 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2326 // Read info about bad chips
2327 fNChips = segm->GetMaximumChipIndex()+1;
2328 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2329 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2330 fChipIsBad = new Bool_t[fNChips];
2331 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2332 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2333 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2334 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2339 //------------------------------------------------------------------------
2340 Double_t AliITStrackerMI::GetEffectiveThickness()
2342 //--------------------------------------------------------------------
2343 // Returns the thickness between the current layer and the vertex (units X0)
2344 //--------------------------------------------------------------------
2347 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2348 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2349 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2353 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2354 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2358 Double_t xn=fgLayers[fI].GetR();
2359 for (Int_t i=0; i<fI; i++) {
2360 Double_t xi=fgLayers[i].GetR();
2361 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2367 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2368 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2371 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2372 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2377 //------------------------------------------------------------------------
2378 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2380 //--------------------------------------------------------------------
2381 // Returns the array of layers between the current layer and the vertex
2382 //--------------------------------------------------------------------
2385 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2386 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2387 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2392 for (int il=fI;il--;) {
2395 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2396 xMS[nl++] = AliITSRecoParam::GetrInsideShield(1);
2399 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2400 xMS[nl++] = AliITSRecoParam::GetrInsideShield(0);
2403 x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2404 xMS[nl++] = fgLayers[il].GetR();
2409 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2410 xMS[nl++] = AliITSRecoParam::GetrPipe();
2416 //------------------------------------------------------------------------
2417 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2418 //-------------------------------------------------------------------
2419 // This function returns number of clusters within the "window"
2420 //--------------------------------------------------------------------
2422 for (Int_t i=fI; i<fN; i++) {
2423 const AliITSRecPoint *c=fClusters[i];
2424 if (c->GetZ() > fZmax) break;
2425 if (c->IsUsed()) continue;
2426 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2427 Double_t y=fR*det.GetPhi() + c->GetY();
2429 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2430 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2432 if (y<fYmin) continue;
2433 if (y>fYmax) continue;
2438 //------------------------------------------------------------------------
2439 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2440 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2442 //--------------------------------------------------------------------
2443 // This function refits the track "track" at the position "x" using
2444 // the clusters from "clusters"
2445 // If "extra"==kTRUE,
2446 // the clusters from overlapped modules get attached to "track"
2447 // If "planeff"==kTRUE,
2448 // special approach for plane efficiency evaluation is applyed
2449 //--------------------------------------------------------------------
2451 Int_t index[AliITSgeomTGeo::kNLayers];
2453 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2454 Int_t nc=clusters->GetNumberOfClusters();
2455 for (k=0; k<nc; k++) {
2456 Int_t idx=clusters->GetClusterIndex(k);
2457 Int_t ilayer=(idx&0xf0000000)>>28;
2461 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2463 //------------------------------------------------------------------------
2464 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2465 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2467 //--------------------------------------------------------------------
2468 // This function refits the track "track" at the position "x" using
2469 // the clusters from array
2470 // If "extra"==kTRUE,
2471 // the clusters from overlapped modules get attached to "track"
2472 // If "planeff"==kTRUE,
2473 // special approach for plane efficiency evaluation is applyed
2474 //--------------------------------------------------------------------
2475 Int_t index[AliITSgeomTGeo::kNLayers];
2477 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2479 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2480 index[k]=clusters[k];
2483 // special for cosmics and TPC prolonged tracks:
2484 // propagate to the innermost of:
2485 // - innermost layer crossed by the track
2486 // - innermost layer where a cluster was associated to the track
2487 static AliITSRecoParam *repa = NULL;
2489 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2491 repa = AliITSRecoParam::GetHighFluxParam();
2492 AliWarning("Using default AliITSRecoParam class");
2495 Int_t evsp=repa->GetEventSpecie();
2497 if(track->GetESDtrack()) trStatus=track->GetStatus();
2498 Int_t innermostlayer=0;
2499 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2501 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2502 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2503 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2504 index[innermostlayer] >= 0 ) break;
2507 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2510 Int_t modstatus=1; // found
2512 Int_t from, to, step;
2513 if (xx > track->GetX()) {
2514 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2517 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2520 TString dir = (step>0 ? "outward" : "inward");
2522 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2523 AliITSlayer &layer=fgLayers[ilayer];
2524 Double_t r=layer.GetR();
2526 if (step<0 && xx>r) break;
2528 // material between SSD and SDD, SDD and SPD
2529 Double_t hI=ilayer-0.5*step;
2530 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2531 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2532 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2533 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2536 Double_t oldGlobXYZ[3];
2537 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2539 // continue if we are already beyond this layer
2540 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2541 if(step>0 && oldGlobR > r) continue; // going outward
2542 if(step<0 && oldGlobR < r) continue; // going inward
2545 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2547 Int_t idet=layer.FindDetectorIndex(phi,z);
2549 // check if we allow a prolongation without point for large-eta tracks
2550 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2552 modstatus = 4; // out in z
2553 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2554 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2557 // apply correction for material of the current layer
2558 // add time if going outward
2559 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2563 if (idet<0) return kFALSE;
2566 const AliITSdetector &det=layer.GetDetector(idet);
2567 // only for ITS-SA tracks refit
2568 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2570 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2572 track->SetDetectorIndex(idet);
2573 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2575 Double_t dz,zmin,zmax,dy,ymin,ymax;
2577 const AliITSRecPoint *clAcc=0;
2578 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2580 Int_t idx=index[ilayer];
2581 if (idx>=0) { // cluster in this layer
2582 modstatus = 6; // no refit
2583 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2585 if (idet != cl->GetDetectorIndex()) {
2586 idet=cl->GetDetectorIndex();
2587 const AliITSdetector &detc=layer.GetDetector(idet);
2588 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2589 track->SetDetectorIndex(idet);
2590 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2592 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2593 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2597 modstatus = 1; // found
2602 } else { // no cluster in this layer
2604 modstatus = 3; // skipped
2605 // Plane Eff determination:
2606 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2607 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2608 UseTrackForPlaneEff(track,ilayer);
2611 modstatus = 5; // no cls in road
2613 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2614 dz = 0.5*(zmax-zmin);
2615 dy = 0.5*(ymax-ymin);
2616 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2617 if (dead==1) modstatus = 7; // holes in z in SPD
2618 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2623 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2624 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2626 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2629 if (extra && clAcc) { // search for extra clusters in overlapped modules
2630 AliITStrackV2 tmp(*track);
2631 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2632 layer.SelectClusters(zmin,zmax,ymin,ymax);
2634 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2636 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2637 Double_t tolerance=0.1;
2638 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2639 // only clusters in another module! (overlaps)
2640 idetExtra = clExtra->GetDetectorIndex();
2641 if (idet == idetExtra) continue;
2643 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2645 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2646 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2647 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2648 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2650 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2651 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2654 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2655 track->SetExtraModule(ilayer,idetExtra);
2657 } // end search for extra clusters in overlapped modules
2659 // Correct for material of the current layer
2661 // add time if going outward
2662 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2663 track->SetCheckInvariant(kTRUE);
2664 } // end loop on layers
2666 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2670 //------------------------------------------------------------------------
2671 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2674 // calculate normalized chi2
2675 // return NormalizedChi2(track,0);
2678 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2679 // track->fdEdxMismatch=0;
2680 Float_t dedxmismatch =0;
2681 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2683 for (Int_t i = 0;i<6;i++){
2684 if (track->GetClIndex(i)>=0){
2685 Float_t cerry, cerrz;
2686 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2688 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2691 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2692 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2693 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2695 cchi2+=(0.5-ratio)*10.;
2696 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2697 dedxmismatch+=(0.5-ratio)*10.;
2701 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2702 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2703 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2704 if (i<2) chi2+=2*cl->GetDeltaProbability();
2710 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2711 track->SetdEdxMismatch(dedxmismatch);
2715 for (Int_t i = 0;i<4;i++){
2716 if (track->GetClIndex(i)>=0){
2717 Float_t cerry, cerrz;
2718 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2719 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2722 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2723 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2727 for (Int_t i = 4;i<6;i++){
2728 if (track->GetClIndex(i)>=0){
2729 Float_t cerry, cerrz;
2730 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2731 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2734 Float_t cerryb, cerrzb;
2735 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2736 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2739 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2740 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2745 if (track->GetESDtrack()->GetTPCsignal()>85){
2746 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2748 chi2+=(0.5-ratio)*5.;
2751 chi2+=(ratio-2.0)*3;
2755 Double_t match = TMath::Sqrt(track->GetChi22());
2756 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2757 if (!track->GetConstrain()) {
2758 if (track->GetNumberOfClusters()>2) {
2759 match/=track->GetNumberOfClusters()-2.;
2764 if (match<0) match=0;
2766 // penalty factor for missing points (NDeadZone>0), but no penalty
2767 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2768 // or there is a dead from OCDB)
2769 Float_t deadzonefactor = 0.;
2770 if (track->GetNDeadZone()>0.) {
2771 Int_t sumDeadZoneProbability=0;
2772 for(Int_t ilay=0;ilay<6;ilay++) {
2773 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2775 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2776 if(nDeadZoneWithProbNot1>0) {
2777 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2778 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2779 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2781 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2782 deadzonefactor = 3.*(1.1-deadZoneProbability);
2786 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2787 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2788 1./(1.+track->GetNSkipped()));
2789 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2790 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2793 //------------------------------------------------------------------------
2794 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2797 // return matching chi2 between two tracks
2798 Double_t largeChi2=1000.;
2800 AliITStrackMI track3(*track2);
2801 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2803 vec(0,0)=track1->GetY() - track3.GetY();
2804 vec(1,0)=track1->GetZ() - track3.GetZ();
2805 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2806 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2807 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2810 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2811 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2812 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2813 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2814 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2816 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2817 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2818 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2819 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2821 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2822 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2823 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2825 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2826 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2828 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2831 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2832 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2835 //------------------------------------------------------------------------
2836 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2839 // return probability that given point (characterized by z position and error)
2840 // is in SPD dead zone
2841 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2843 Double_t probability = 0.;
2844 Double_t nearestz = 0.,distz=0.;
2845 Int_t nearestzone = -1;
2846 Double_t mindistz = 1000.;
2848 // find closest dead zone
2849 for (Int_t i=0; i<3; i++) {
2850 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2851 if (distz<mindistz) {
2853 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2858 // too far from dead zone
2859 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2862 Double_t zmin, zmax;
2863 if (nearestzone==0) { // dead zone at z = -7
2864 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2865 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2866 } else if (nearestzone==1) { // dead zone at z = 0
2867 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2868 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2869 } else if (nearestzone==2) { // dead zone at z = +7
2870 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2871 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2876 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2878 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2879 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2880 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2883 //------------------------------------------------------------------------
2884 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2887 // calculate normalized chi2
2889 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2891 for (Int_t i = 0;i<6;i++){
2892 if (TMath::Abs(track->GetDy(i))>0){
2893 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2894 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2897 else{chi2[i]=10000;}
2900 TMath::Sort(6,chi2,index,kFALSE);
2901 Float_t max = float(ncl)*fac-1.;
2902 Float_t sumchi=0, sumweight=0;
2903 for (Int_t i=0;i<max+1;i++){
2904 Float_t weight = (i<max)?1.:(max+1.-i);
2905 sumchi+=weight*chi2[index[i]];
2908 Double_t normchi2 = sumchi/sumweight;
2911 //------------------------------------------------------------------------
2912 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2915 // calculate normalized chi2
2916 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2919 for (Int_t i=0;i<6;i++){
2920 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2921 Double_t sy1 = forwardtrack->GetSigmaY(i);
2922 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2923 Double_t sy2 = backtrack->GetSigmaY(i);
2924 Double_t sz2 = backtrack->GetSigmaZ(i);
2925 if (i<2){ sy2=1000.;sz2=1000;}
2927 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2928 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2930 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2931 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2933 res+= nz0*nz0+ny0*ny0;
2936 if (npoints>1) return
2937 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2938 //2*forwardtrack->fNUsed+
2939 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2940 1./(1.+forwardtrack->GetNSkipped()));
2943 //------------------------------------------------------------------------
2944 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2945 //--------------------------------------------------------------------
2946 // Return pointer to a given cluster
2947 //--------------------------------------------------------------------
2948 Int_t l=(index & 0xf0000000) >> 28;
2949 Int_t c=(index & 0x0fffffff) >> 00;
2950 return fgLayers[l].GetWeight(c);
2952 //------------------------------------------------------------------------
2953 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2955 //---------------------------------------------
2956 // register track to the list
2958 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2961 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2962 Int_t index = track->GetClusterIndex(icluster);
2963 Int_t l=(index & 0xf0000000) >> 28;
2964 Int_t c=(index & 0x0fffffff) >> 00;
2965 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2966 for (Int_t itrack=0;itrack<4;itrack++){
2967 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2968 fgLayers[l].SetClusterTracks(itrack,c,id);
2974 //------------------------------------------------------------------------
2975 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2977 //---------------------------------------------
2978 // unregister track from the list
2979 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2980 Int_t index = track->GetClusterIndex(icluster);
2981 Int_t l=(index & 0xf0000000) >> 28;
2982 Int_t c=(index & 0x0fffffff) >> 00;
2983 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2984 for (Int_t itrack=0;itrack<4;itrack++){
2985 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2986 fgLayers[l].SetClusterTracks(itrack,c,-1);
2991 //------------------------------------------------------------------------
2992 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2994 //-------------------------------------------------------------
2995 //get number of shared clusters
2996 //-------------------------------------------------------------
2998 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2999 // mean number of clusters
3000 Float_t *ny = GetNy(id), *nz = GetNz(id);
3003 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
3004 Int_t index = track->GetClusterIndex(icluster);
3005 Int_t l=(index & 0xf0000000) >> 28;
3006 Int_t c=(index & 0x0fffffff) >> 00;
3007 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3008 // if (ny[l]<1.e-13){
3009 // printf("problem\n");
3011 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3015 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3016 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3017 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3019 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3022 deltan = (cl->GetNz()-nz[l]);
3024 if (deltan>2.0) continue; // extended - highly probable shared cluster
3025 weight = 2./TMath::Max(3.+deltan,2.);
3027 for (Int_t itrack=0;itrack<4;itrack++){
3028 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3030 clist[l] = (AliITSRecPoint*)GetCluster(index);
3031 track->SetSharedWeight(l,weight);
3037 track->SetNUsed(shared);
3040 //------------------------------------------------------------------------
3041 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3044 // find first shared track
3046 // mean number of clusters
3047 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
3049 for (Int_t i=0;i<6;i++) overlist[i]=-1;
3050 Int_t sharedtrack=100000;
3051 Int_t tracks[24],trackindex=0;
3052 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3054 for (Int_t icluster=0;icluster<6;icluster++){
3055 if (clusterlist[icluster]<0) continue;
3056 Int_t index = clusterlist[icluster];
3057 Int_t l=(index & 0xf0000000) >> 28;
3058 Int_t c=(index & 0x0fffffff) >> 00;
3059 // if (ny[l]<1.e-13){
3060 // printf("problem\n");
3062 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3063 //if (l>3) continue;
3064 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3067 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3068 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3069 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3071 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3074 deltan = (cl->GetNz()-nz[l]);
3076 if (deltan>2.0) continue; // extended - highly probable shared cluster
3078 for (Int_t itrack=3;itrack>=0;itrack--){
3079 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3080 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3081 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
3086 if (trackindex==0) return -1;
3088 sharedtrack = tracks[0];
3090 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3093 Int_t tracks2[24], cluster[24];
3094 for (Int_t i=0;i<24;i++){ tracks2[i]=-1; cluster[i]=0;}
3097 for (Int_t i=0;i<trackindex;i++){
3098 if (tracks[i]<0) continue;
3099 tracks2[index] = tracks[i];
3101 for (Int_t j=i+1;j<trackindex;j++){
3102 if (tracks[j]<0) continue;
3103 if (tracks[j]==tracks[i]){
3111 for (Int_t i=0;i<index;i++){
3112 if (cluster[index]>max) {
3113 sharedtrack=tracks2[index];
3120 if (sharedtrack>=100000) return -1;
3122 // make list of overlaps
3124 for (Int_t icluster=0;icluster<6;icluster++){
3125 if (clusterlist[icluster]<0) continue;
3126 Int_t index = clusterlist[icluster];
3127 Int_t l=(index & 0xf0000000) >> 28;
3128 Int_t c=(index & 0x0fffffff) >> 00;
3129 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3130 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3132 if (cl->GetNy()>2) continue;
3133 if (cl->GetNz()>2) continue;
3136 if (cl->GetNy()>3) continue;
3137 if (cl->GetNz()>3) continue;
3140 for (Int_t itrack=3;itrack>=0;itrack--){
3141 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3142 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3150 //------------------------------------------------------------------------
3151 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3153 // try to find track hypothesys without conflicts
3154 // with minimal chi2;
3155 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3156 Int_t entries1 = arr1->GetEntriesFast();
3157 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3158 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3159 Int_t entries2 = arr2->GetEntriesFast();
3160 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3162 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3163 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3164 if (track10->Pt()>0.5+track20->Pt()) return track10;
3166 for (Int_t itrack=0;itrack<entries1;itrack++){
3167 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3168 UnRegisterClusterTracks(track,trackID1);
3171 for (Int_t itrack=0;itrack<entries2;itrack++){
3172 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3173 UnRegisterClusterTracks(track,trackID2);
3177 Float_t maxconflicts=6;
3178 Double_t maxchi2 =1000.;
3180 // get weight of hypothesys - using best hypothesy
3183 Int_t list1[6],list2[6];
3184 AliITSRecPoint *clist1[6], *clist2[6] ;
3185 RegisterClusterTracks(track10,trackID1);
3186 RegisterClusterTracks(track20,trackID2);
3187 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3188 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3189 UnRegisterClusterTracks(track10,trackID1);
3190 UnRegisterClusterTracks(track20,trackID2);
3193 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3194 Float_t nerry[6],nerrz[6];
3195 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3196 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3197 for (Int_t i=0;i<6;i++){
3198 if ( (erry1[i]>0) && (erry2[i]>0)) {
3199 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3200 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3202 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3203 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3205 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3206 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3207 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3210 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3211 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3212 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3220 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3221 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3222 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3223 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3225 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3226 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3227 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3229 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3230 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3231 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3234 Double_t sumw = w1+w2;
3238 w1 = (d2+0.5)/(d1+d2+1.);
3239 w2 = (d1+0.5)/(d1+d2+1.);
3241 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3242 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3244 // get pair of "best" hypothesys
3246 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3247 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3249 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3250 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3251 //if (track1->fFakeRatio>0) continue;
3252 RegisterClusterTracks(track1,trackID1);
3253 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3254 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3256 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3257 //if (track2->fFakeRatio>0) continue;
3259 RegisterClusterTracks(track2,trackID2);
3260 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3261 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3262 UnRegisterClusterTracks(track2,trackID2);
3264 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3265 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3266 if (nskipped>0.5) continue;
3268 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3269 if (conflict1+1<cconflict1) continue;
3270 if (conflict2+1<cconflict2) continue;
3274 for (Int_t i=0;i<6;i++){
3276 Float_t c1 =0.; // conflict coeficients
3278 if (clist1[i]&&clist2[i]){
3281 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3284 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3286 c1 = 2./TMath::Max(3.+deltan,2.);
3287 c2 = 2./TMath::Max(3.+deltan,2.);
3293 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3296 deltan = (clist1[i]->GetNz()-nz1[i]);
3298 c1 = 2./TMath::Max(3.+deltan,2.);
3305 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3308 deltan = (clist2[i]->GetNz()-nz2[i]);
3310 c2 = 2./TMath::Max(3.+deltan,2.);
3316 if (TMath::Abs(track1->GetDy(i))>0.) {
3317 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3318 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3319 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3320 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3322 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3325 if (TMath::Abs(track2->GetDy(i))>0.) {
3326 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3327 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3328 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3329 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3332 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3334 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3335 if (chi21>0) sum+=w1;
3336 if (chi22>0) sum+=w2;
3339 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3340 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3341 Double_t normchi2 = 2*conflict+sumchi2/sum;
3342 if ( normchi2 <maxchi2 ){
3345 maxconflicts = conflict;
3349 UnRegisterClusterTracks(track1,trackID1);
3352 // if (maxconflicts<4 && maxchi2<th0){
3353 if (maxchi2<th0*2.){
3354 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3355 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3356 track1->SetChi2MIP(5,maxconflicts);
3357 track1->SetChi2MIP(6,maxchi2);
3358 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3359 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3360 track1->SetChi2MIP(8,index1);
3361 fBestTrackIndex[trackID1] =index1;
3362 UpdateESDtrack(track1, AliESDtrack::kITSin);
3363 original->SetWinner(track1);
3365 else if (track10->GetChi2MIP(0)<th1){
3366 track10->SetChi2MIP(5,maxconflicts);
3367 track10->SetChi2MIP(6,maxchi2);
3368 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3369 UpdateESDtrack(track10,AliESDtrack::kITSin);
3370 original->SetWinner(track10);
3373 for (Int_t itrack=0;itrack<entries1;itrack++){
3374 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3375 UnRegisterClusterTracks(track,trackID1);
3378 for (Int_t itrack=0;itrack<entries2;itrack++){
3379 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3380 UnRegisterClusterTracks(track,trackID2);
3383 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3384 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3385 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3386 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3387 RegisterClusterTracks(track10,trackID1);
3389 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3390 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3391 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3392 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3393 RegisterClusterTracks(track20,trackID2);
3398 //------------------------------------------------------------------------
3399 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3400 //--------------------------------------------------------------------
3401 // This function marks clusters assigned to the track
3402 //--------------------------------------------------------------------
3403 AliTracker::UseClusters(t,from);
3405 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3406 //if (c->GetQ()>2) c->Use();
3407 if (c->GetSigmaZ2()>0.1) c->Use();
3408 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3409 //if (c->GetQ()>2) c->Use();
3410 if (c->GetSigmaZ2()>0.1) c->Use();
3413 //------------------------------------------------------------------------
3414 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3416 //------------------------------------------------------------------
3417 // add track to the list of hypothesys
3418 //------------------------------------------------------------------
3421 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3423 array = new TObjArray(10);
3424 fTrackHypothesys.AddAtAndExpand(array,esdindex);
3426 array->AddLast(track);
3428 //------------------------------------------------------------------------
3429 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3431 //-------------------------------------------------------------------
3432 // compress array of track hypothesys
3433 // keep only maxsize best hypothesys
3434 //-------------------------------------------------------------------
3435 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3436 if (! (fTrackHypothesys.At(esdindex)) ) return;
3437 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3438 Int_t entries = array->GetEntriesFast();
3440 //- find preliminary besttrack as a reference
3441 Float_t minchi2=10000;
3443 AliITStrackMI * besttrack=0;
3445 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3446 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3447 if (!track) continue;
3448 Float_t chi2 = NormalizedChi2(track,0);
3450 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3451 track->SetLabel(tpcLabel);
3453 track->SetFakeRatio(1.);
3454 CookLabel(track,0.); //For comparison only
3456 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3457 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3458 if (track->GetNumberOfClusters()<maxn) continue;
3459 maxn = track->GetNumberOfClusters();
3460 // if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3467 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3468 delete array->RemoveAt(itrack);
3472 if (!besttrack) return;
3475 //take errors of best track as a reference
3476 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3477 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3478 for (Int_t j=0;j<6;j++) {
3479 if (besttrack->GetClIndex(j)>=0){
3480 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3481 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3482 ny[j] = besttrack->GetNy(j);
3483 nz[j] = besttrack->GetNz(j);
3487 // calculate normalized chi2
3489 Float_t * chi2 = new Float_t[entries];
3490 Int_t * index = new Int_t[entries];
3491 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3492 for (Int_t itrack=0;itrack<entries;itrack++){
3493 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3495 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3496 double chi2t = GetNormalizedChi2(track, mode);
3497 track->SetChi2MIP(0,chi2t);
3498 if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3499 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3500 chi2[itrack] = chi2t;
3503 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3504 delete array->RemoveAt(itrack);
3510 TMath::Sort(entries,chi2,index,kFALSE);
3511 besttrack = (AliITStrackMI*)array->At(index[0]);
3512 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3513 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3514 for (Int_t j=0;j<6;j++){
3515 if (besttrack->GetClIndex(j)>=0){
3516 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3517 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3518 ny[j] = besttrack->GetNy(j);
3519 nz[j] = besttrack->GetNz(j);
3524 // calculate one more time with updated normalized errors
3525 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3526 for (Int_t itrack=0;itrack<entries;itrack++){
3527 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3529 double chi2t = GetNormalizedChi2(track, mode);
3530 track->SetChi2MIP(0,chi2t);
3531 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3532 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3533 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3534 chi2[itrack] = chi2t; //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3537 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3538 delete array->RemoveAt(itrack);
3543 entries = array->GetEntriesFast();
3547 TObjArray * newarray = new TObjArray();
3548 TMath::Sort(entries,chi2,index,kFALSE);
3549 besttrack = (AliITStrackMI*)array->At(index[0]);
3551 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3553 for (Int_t j=0;j<6;j++){
3554 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3555 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3556 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3557 ny[j] = besttrack->GetNy(j);
3558 nz[j] = besttrack->GetNz(j);
3561 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3562 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3563 Float_t minn = besttrack->GetNumberOfClusters()-3;
3565 for (Int_t i=0;i<entries;i++){
3566 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3567 if (!track) continue;
3568 if (accepted>maxcut) break;
3569 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3570 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3571 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3572 delete array->RemoveAt(index[i]);
3576 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3577 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3578 if (!shortbest) accepted++;
3580 newarray->AddLast(array->RemoveAt(index[i]));
3581 for (Int_t j=0;j<6;j++){
3583 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3584 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3585 ny[j] = track->GetNy(j);
3586 nz[j] = track->GetNz(j);
3591 delete array->RemoveAt(index[i]);
3595 delete fTrackHypothesys.RemoveAt(esdindex);
3596 fTrackHypothesys.AddAt(newarray,esdindex);
3600 delete fTrackHypothesys.RemoveAt(esdindex);
3606 //------------------------------------------------------------------------
3607 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3609 //-------------------------------------------------------------
3610 // try to find best hypothesy
3611 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3612 // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3613 //-------------------------------------------------------------
3614 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3615 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3616 if (!array) return 0;
3617 Int_t entries = array->GetEntriesFast();
3618 if (!entries) return 0;
3619 Float_t minchi2 = 100000;
3620 AliITStrackMI * besttrack=0;
3622 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3623 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3624 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3625 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3627 for (Int_t i=0;i<entries;i++){
3628 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3629 if (!track) continue;
3630 Float_t sigmarfi,sigmaz;
3631 GetDCASigma(track,sigmarfi,sigmaz);
3632 track->SetDnorm(0,sigmarfi);
3633 track->SetDnorm(1,sigmaz);
3635 track->SetChi2MIP(1,1000000);
3636 track->SetChi2MIP(2,1000000);
3637 track->SetChi2MIP(3,1000000);
3640 backtrack = new(backtrack) AliITStrackMI(*track);
3641 if (track->GetConstrain()) {
3642 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3643 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3644 if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3645 else {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3647 backtrack->ResetCovariance(10.);
3649 backtrack->ResetCovariance(10.);
3651 backtrack->ResetClusters();
3653 Double_t x = original->GetX();
3654 if (!RefitAt(x,backtrack,track)) continue;
3656 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3657 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3658 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3659 track->SetChi22(GetMatchingChi2(backtrack,original));
3661 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3662 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3663 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3666 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3668 //forward track - without constraint
3669 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3670 forwardtrack->ResetClusters();
3672 if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue; // w/o fwd track MIP03 is meaningless
3673 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3674 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3675 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3677 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3678 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3679 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3680 forwardtrack->SetD(0,track->GetD(0));
3681 forwardtrack->SetD(1,track->GetD(1));
3684 AliITSRecPoint* clist[6];
3685 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3686 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3689 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3690 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3691 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3692 track->SetChi2MIP(3,1000);
3695 Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed(); //RS
3696 if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3697 else chi2 += track->GetNUsed();
3699 for (Int_t ichi=0;ichi<5;ichi++){
3700 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3702 if (chi2 < minchi2){
3703 //besttrack = new AliITStrackMI(*forwardtrack);
3705 besttrack->SetLabel(track->GetLabel());
3706 besttrack->SetFakeRatio(track->GetFakeRatio());
3708 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3709 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3710 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3714 delete forwardtrack;
3716 if (!besttrack) return 0;
3719 for (Int_t i=0;i<entries;i++){
3720 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3722 if (!track) continue;
3724 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3725 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3726 // RS: don't apply this cut when fSelectBestMIP03 is on
3727 || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3729 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3730 delete array->RemoveAt(i);
3740 SortTrackHypothesys(esdindex,checkmax,1);
3742 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3743 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3744 besttrack = (AliITStrackMI*)array->At(0);
3745 if (!besttrack) return 0;
3746 besttrack->SetChi2MIP(8,0);
3747 fBestTrackIndex[esdindex]=0;
3748 entries = array->GetEntriesFast();
3749 AliITStrackMI *longtrack =0;
3751 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3752 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3753 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3754 if (!track->GetConstrain()) continue;
3755 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3756 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3757 if (track->GetChi2MIP(0)>4.) continue;
3758 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3761 //if (longtrack) besttrack=longtrack;
3763 // RS do shared cluster analysis here only if the new sharing analysis is not requested
3764 //RRR if (fFlagFakes) return besttrack;
3767 AliITSRecPoint * clist[6];
3768 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3769 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3770 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3771 RegisterClusterTracks(besttrack,esdindex);
3778 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3779 if (sharedtrack>=0){
3781 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3783 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3789 if (shared>2.5) return 0;
3790 if (shared>1.0) return besttrack;
3792 // Don't sign clusters if not gold track
3794 if (!besttrack->IsGoldPrimary()) return besttrack;
3795 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3797 if (fConstraint[fPass]){
3801 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3802 for (Int_t i=0;i<6;i++){
3803 Int_t index = besttrack->GetClIndex(i);
3804 if (index<0) continue;
3805 Int_t ilayer = (index & 0xf0000000) >> 28;
3806 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3807 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3809 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3810 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3811 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3812 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3813 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3814 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3816 Bool_t cansign = kTRUE;
3817 for (Int_t itrack=0;itrack<entries; itrack++){
3818 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3819 if (!track) continue;
3820 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3821 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3827 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3828 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3829 if (!c->IsUsed()) c->Use();
3835 //------------------------------------------------------------------------
3836 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3839 // get "best" hypothesys
3842 Int_t nentries = itsTracks.GetEntriesFast();
3843 for (Int_t i=0;i<nentries;i++){
3844 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3845 if (!track) continue;
3846 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3847 if (!array) continue;
3848 if (array->GetEntriesFast()<=0) continue;
3850 AliITStrackMI* longtrack=0;
3852 Float_t maxchi2=1000;
3853 for (Int_t j=0;j<array->GetEntriesFast();j++){
3854 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3855 if (!trackHyp) continue;
3856 if (trackHyp->GetGoldV0()) {
3857 longtrack = trackHyp; //gold V0 track taken
3860 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3861 Float_t chi2 = trackHyp->GetChi2MIP(0);
3862 if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3863 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3865 if (fAfterV0){ // ??? RS
3866 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3868 if (chi2 > maxchi2) continue;
3869 minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3870 if (fSelectBestMIP03) minn++; // allow next to longest to win
3877 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3878 if (!longtrack) {longtrack = besttrack;}
3879 else besttrack= longtrack;
3883 AliITSRecPoint * clist[6];
3884 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3886 track->SetNUsed(shared);
3887 track->SetNSkipped(besttrack->GetNSkipped());
3888 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3890 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3894 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3895 //if (sharedtrack==-1) sharedtrack=0;
3896 if (sharedtrack>=0) {
3897 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3900 if (besttrack&&fAfterV0) {
3901 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3902 track->SetWinner(besttrack);
3905 if (fConstraint[fPass]) {
3906 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3907 track->SetWinner(besttrack);
3909 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3910 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3911 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3918 //------------------------------------------------------------------------
3919 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3922 // RS: flag those tracks which are suxpected to have fake clusters
3924 const double kThreshPt = 0.5;
3925 AliRefArray *refArr[6];
3927 for (int i=0;i<6;i++) {
3928 int ncl = fgLayers[i].GetNumberOfClusters();
3929 refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3931 Int_t nentries = itsTracks.GetEntriesFast();
3933 // fill cluster->track associations
3934 for (Int_t itr=0;itr<nentries;itr++){
3935 AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3936 if (!track) continue;
3937 AliITStrackMI* trackITS = track->GetWinner();
3938 if (!trackITS) continue;
3939 for (int il=trackITS->GetNumberOfClusters();il--;) {
3940 int idx = trackITS->GetClusterIndex(il);
3941 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3942 // if (c>fgLayers[l].GetNumberOfClusters()) continue;
3943 refArr[l]->AddReference(c, itr);
3947 const UInt_t kMaxRef = 100;
3948 UInt_t crefs[kMaxRef];
3950 // process tracks with shared clusters
3951 for (int itr=0;itr<nentries;itr++){
3952 AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3953 AliITStrackMI* trackH0 = track0->GetWinner();
3954 if (!trackH0) continue;
3955 AliESDtrack* esd0 = track0->GetESDtrack();
3957 for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3958 int idx = trackH0->GetClusterIndex(il);
3959 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3960 ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3961 if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3962 esd0->SetITSSharedFlag(l);
3963 for (int ir=ncrefs;ir--;) {
3964 if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3965 AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3966 AliITStrackMI* trackH1 = track1->GetWinner();
3967 AliESDtrack* esd1 = track1->GetESDtrack();
3968 esd1->SetITSSharedFlag(l);
3970 double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3971 if (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3972 else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3974 // select the one with smallest chi2's product
3975 res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3976 res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3978 if (res<0) esd1->SetITSFakeFlag(); // esd0 is winner
3979 else esd0->SetITSFakeFlag(); // esd1 is winner
3986 for (int i=6;i--;) delete refArr[i];
3989 //------------------------------------------------------------------------
3990 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3991 //--------------------------------------------------------------------
3992 //This function "cooks" a track label. If label<0, this track is fake.
3993 //--------------------------------------------------------------------
3996 if (track->GetESDtrack()){
3997 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3998 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3999 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
4001 track->SetChi2MIP(9,0);
4003 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4004 Int_t cindex = track->GetClusterIndex(i);
4005 Int_t l=(cindex & 0xf0000000) >> 28;
4006 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4008 for (Int_t ind=0;ind<3;ind++){
4009 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
4010 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
4012 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4015 Int_t nclusters = track->GetNumberOfClusters();
4016 if (nclusters > 0) //PH Some tracks don't have any cluster
4017 track->SetFakeRatio(double(nwrong)/double(nclusters));
4018 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4019 track->SetLabel(-tpcLabel);
4021 track->SetLabel(tpcLabel);
4023 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4026 //------------------------------------------------------------------------
4027 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4029 // Fill the dE/dx in this track
4031 track->SetChi2MIP(9,0);
4032 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4033 Int_t cindex = track->GetClusterIndex(i);
4034 Int_t l=(cindex & 0xf0000000) >> 28;
4035 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4036 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4038 for (Int_t ind=0;ind<3;ind++){
4039 if (cl->GetLabel(ind)==lab) isWrong=0;
4041 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4045 track->CookdEdx(low,up);
4047 //------------------------------------------------------------------------
4048 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4050 // Create some arrays
4052 if (fCoefficients) delete []fCoefficients;
4053 fCoefficients = new Float_t[ntracks*48];
4054 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4056 //------------------------------------------------------------------------
4057 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4060 // Compute predicted chi2
4062 // Take into account the mis-alignment (bring track to cluster plane)
4063 Double_t xTrOrig=track->GetX();
4064 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4065 Float_t erry,errz,covyz;
4066 Float_t theta = track->GetTgl();
4067 Float_t phi = track->GetSnp();
4068 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4069 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4070 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()));
4071 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()));
4072 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4073 // Bring the track back to detector plane in ideal geometry
4074 // [mis-alignment will be accounted for in UpdateMI()]
4075 if (!track->Propagate(xTrOrig)) return 1000.;
4077 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4078 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4080 chi2+=0.5*TMath::Min(delta/2,2.);
4081 chi2+=2.*cluster->GetDeltaProbability();
4084 track->SetNy(layer,ny);
4085 track->SetNz(layer,nz);
4086 track->SetSigmaY(layer,erry);
4087 track->SetSigmaZ(layer, errz);
4088 track->SetSigmaYZ(layer,covyz);
4089 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4090 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4094 //------------------------------------------------------------------------
4095 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4100 Int_t layer = (index & 0xf0000000) >> 28;
4101 track->SetClIndex(layer, index);
4102 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4103 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4104 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4105 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4109 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
4112 // Take into account the mis-alignment (bring track to cluster plane)
4113 Double_t xTrOrig=track->GetX();
4114 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4115 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4116 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4117 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
4119 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4122 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4123 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4124 c.SetSigmaYZ(track->GetSigmaYZ(layer));
4127 Int_t updated = track->UpdateMI(&c,chi2,index);
4128 // Bring the track back to detector plane in ideal geometry
4129 if (!track->Propagate(xTrOrig)) return 0;
4131 if(!updated) AliDebug(2,"update failed");
4135 //------------------------------------------------------------------------
4136 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4139 //DCA sigmas parameterization
4140 //to be paramterized using external parameters in future
4143 Double_t curv=track->GetC();
4144 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4145 sigmaz = 0.0110+4.37*TMath::Abs(curv);
4147 //------------------------------------------------------------------------
4148 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4151 // Clusters from delta electrons?
4153 Int_t entries = clusterArray->GetEntriesFast();
4154 if (entries<4) return;
4155 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4156 Int_t layer = cluster->GetLayer();
4157 if (layer>1) return;
4159 Int_t ncandidates=0;
4160 Float_t r = (layer>0)? 7:4;
4162 for (Int_t i=0;i<entries;i++){
4163 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4164 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4165 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4166 index[ncandidates] = i; //candidate to belong to delta electron track
4168 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4169 cl0->SetDeltaProbability(1);
4175 for (Int_t i=0;i<ncandidates;i++){
4176 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4177 if (cl0->GetDeltaProbability()>0.8) continue;
4180 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4181 sumy=sumz=sumy2=sumyz=sumw=0.0;
4182 for (Int_t j=0;j<ncandidates;j++){
4184 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4186 Float_t dz = cl0->GetZ()-cl1->GetZ();
4187 Float_t dy = cl0->GetY()-cl1->GetY();
4188 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4189 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4190 y[ncl] = cl1->GetY();
4191 z[ncl] = cl1->GetZ();
4192 sumy+= y[ncl]*weight;
4193 sumz+= z[ncl]*weight;
4194 sumy2+=y[ncl]*y[ncl]*weight;
4195 sumyz+=y[ncl]*z[ncl]*weight;
4200 if (ncl<4) continue;
4201 Float_t det = sumw*sumy2 - sumy*sumy;
4203 if (TMath::Abs(det)>0.01){
4204 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
4205 Float_t k = (sumyz*sumw - sumy*sumz)/det;
4206 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4209 Float_t z0 = sumyz/sumy;
4210 delta = TMath::Abs(cl0->GetZ()-z0);
4213 cl0->SetDeltaProbability(1-20.*delta);
4217 //------------------------------------------------------------------------
4218 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4223 track->UpdateESDtrack(flags);
4224 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4225 if (oldtrack) delete oldtrack;
4226 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4227 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4228 // printf("Problem\n");
4231 //------------------------------------------------------------------------
4232 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4234 // Get nearest upper layer close to the point xr.
4235 // rough approximation
4237 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4238 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4240 for (Int_t i=0;i<6;i++){
4241 if (radius<kRadiuses[i]){
4248 //------------------------------------------------------------------------
4249 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4250 //--------------------------------------------------------------------
4251 // Fill a look-up table with mean material
4252 //--------------------------------------------------------------------
4256 Double_t point1[3],point2[3];
4257 Double_t phi,cosphi,sinphi,z;
4258 // 0-5 layers, 6 pipe, 7-8 shields
4259 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4260 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4262 Int_t ifirst=0,ilast=0;
4263 if(material.Contains("Pipe")) {
4265 } else if(material.Contains("Shields")) {
4267 } else if(material.Contains("Layers")) {
4270 Error("BuildMaterialLUT","Wrong layer name\n");
4273 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4274 Double_t param[5]={0.,0.,0.,0.,0.};
4275 for (Int_t i=0; i<n; i++) {
4276 phi = 2.*TMath::Pi()*gRandom->Rndm();
4277 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4278 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4279 point1[0] = rmin[imat]*cosphi;
4280 point1[1] = rmin[imat]*sinphi;
4282 point2[0] = rmax[imat]*cosphi;
4283 point2[1] = rmax[imat]*sinphi;
4285 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4286 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4288 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4290 fxOverX0Layer[imat] = param[1];
4291 fxTimesRhoLayer[imat] = param[0]*param[4];
4292 } else if(imat==6) {
4293 fxOverX0Pipe = param[1];
4294 fxTimesRhoPipe = param[0]*param[4];
4295 } else if(imat==7) {
4296 fxOverX0Shield[0] = param[1];
4297 fxTimesRhoShield[0] = param[0]*param[4];
4298 } else if(imat==8) {
4299 fxOverX0Shield[1] = param[1];
4300 fxTimesRhoShield[1] = param[0]*param[4];
4304 printf("%s\n",material.Data());
4305 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4306 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4307 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4308 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4309 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4310 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4311 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4312 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4313 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4317 //------------------------------------------------------------------------
4318 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4319 TString direction) {
4320 //-------------------------------------------------------------------
4321 // Propagate beyond beam pipe and correct for material
4322 // (material budget in different ways according to fUseTGeo value)
4323 // Add time if going outward (PropagateTo or PropagateToTGeo)
4324 //-------------------------------------------------------------------
4326 // Define budget mode:
4327 // 0: material from AliITSRecoParam (hard coded)
4328 // 1: material from TGeo in one step (on the fly)
4329 // 2: material from lut
4330 // 3: material from TGeo in one step (same for all hypotheses)
4343 if(fTrackingPhase.Contains("Clusters2Tracks"))
4344 { mode=3; } else { mode=1; }
4347 if(fTrackingPhase.Contains("Clusters2Tracks"))
4348 { mode=3; } else { mode=2; }
4354 if(fTrackingPhase.Contains("Default")) mode=0;
4356 Int_t index=fCurrentEsdTrack;
4358 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4359 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4361 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4363 Double_t xOverX0,x0,lengthTimesMeanDensity;
4367 xOverX0 = AliITSRecoParam::GetdPipe();
4368 x0 = AliITSRecoParam::GetX0Be();
4369 lengthTimesMeanDensity = xOverX0*x0;
4370 lengthTimesMeanDensity *= dir;
4371 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4374 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4377 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4378 xOverX0 = fxOverX0Pipe;
4379 lengthTimesMeanDensity = fxTimesRhoPipe;
4380 lengthTimesMeanDensity *= dir;
4381 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4384 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4385 if(fxOverX0PipeTrks[index]<0) {
4386 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4387 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4388 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4389 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4390 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4393 xOverX0 = fxOverX0PipeTrks[index];
4394 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4395 lengthTimesMeanDensity *= dir;
4396 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4402 //------------------------------------------------------------------------
4403 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4405 TString direction) {
4406 //-------------------------------------------------------------------
4407 // Propagate beyond SPD or SDD shield and correct for material
4408 // (material budget in different ways according to fUseTGeo value)
4409 // Add time if going outward (PropagateTo or PropagateToTGeo)
4410 //-------------------------------------------------------------------
4412 // Define budget mode:
4413 // 0: material from AliITSRecoParam (hard coded)
4414 // 1: material from TGeo in steps of X cm (on the fly)
4415 // X = AliITSRecoParam::GetStepSizeTGeo()
4416 // 2: material from lut
4417 // 3: material from TGeo in one step (same for all hypotheses)
4430 if(fTrackingPhase.Contains("Clusters2Tracks"))
4431 { mode=3; } else { mode=1; }
4434 if(fTrackingPhase.Contains("Clusters2Tracks"))
4435 { mode=3; } else { mode=2; }
4441 if(fTrackingPhase.Contains("Default")) mode=0;
4443 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4445 Int_t shieldindex=0;
4446 if (shield.Contains("SDD")) { // SDDouter
4447 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4449 } else if (shield.Contains("SPD")) { // SPDouter
4450 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4453 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4457 // do nothing if we are already beyond the shield
4458 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4459 if(dir<0 && rTrack > rToGo) return 1; // going outward
4460 if(dir>0 && rTrack < rToGo) return 1; // going inward
4464 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4466 Int_t index=2*fCurrentEsdTrack+shieldindex;
4468 Double_t xOverX0,x0,lengthTimesMeanDensity;
4473 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4474 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4475 lengthTimesMeanDensity = xOverX0*x0;
4476 lengthTimesMeanDensity *= dir;
4477 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4480 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4481 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4484 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4485 xOverX0 = fxOverX0Shield[shieldindex];
4486 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4487 lengthTimesMeanDensity *= dir;
4488 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4491 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4492 if(fxOverX0ShieldTrks[index]<0) {
4493 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4494 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4495 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4496 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4497 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4500 xOverX0 = fxOverX0ShieldTrks[index];
4501 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4502 lengthTimesMeanDensity *= dir;
4503 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4509 //------------------------------------------------------------------------
4510 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4512 Double_t oldGlobXYZ[3],
4513 TString direction) {
4514 //-------------------------------------------------------------------
4515 // Propagate beyond layer and correct for material
4516 // (material budget in different ways according to fUseTGeo value)
4517 // Add time if going outward (PropagateTo or PropagateToTGeo)
4518 //-------------------------------------------------------------------
4520 // Define budget mode:
4521 // 0: material from AliITSRecoParam (hard coded)
4522 // 1: material from TGeo in stepsof X cm (on the fly)
4523 // X = AliITSRecoParam::GetStepSizeTGeo()
4524 // 2: material from lut
4525 // 3: material from TGeo in one step (same for all hypotheses)
4538 if(fTrackingPhase.Contains("Clusters2Tracks"))
4539 { mode=3; } else { mode=1; }
4542 if(fTrackingPhase.Contains("Clusters2Tracks"))
4543 { mode=3; } else { mode=2; }
4549 if(fTrackingPhase.Contains("Default")) mode=0;
4551 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4553 Double_t r=fgLayers[layerindex].GetR();
4554 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4556 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4558 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4560 Int_t index=6*fCurrentEsdTrack+layerindex;
4563 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4566 // back before material (no correction)
4568 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4569 if (!t->GetLocalXat(rOld,xOld)) return 0;
4570 if (!t->Propagate(xOld)) return 0;
4574 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4575 lengthTimesMeanDensity = xOverX0*x0;
4576 lengthTimesMeanDensity *= dir;
4577 // Bring the track beyond the material
4578 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4581 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4582 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4585 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4586 xOverX0 = fxOverX0Layer[layerindex];
4587 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4588 lengthTimesMeanDensity *= dir;
4589 // Bring the track beyond the material
4590 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4593 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4594 if(fxOverX0LayerTrks[index]<0) {
4595 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4596 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4597 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4598 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4599 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4600 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4603 xOverX0 = fxOverX0LayerTrks[index];
4604 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4605 lengthTimesMeanDensity *= dir;
4606 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4613 //------------------------------------------------------------------------
4614 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4615 //-----------------------------------------------------------------
4616 // Initialize LUT for storing material for each prolonged track
4617 //-----------------------------------------------------------------
4618 fxOverX0PipeTrks = new Float_t[ntracks];
4619 fxTimesRhoPipeTrks = new Float_t[ntracks];
4620 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4621 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4622 fxOverX0LayerTrks = new Float_t[ntracks*6];
4623 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4625 for(Int_t i=0; i<ntracks; i++) {
4626 fxOverX0PipeTrks[i] = -1.;
4627 fxTimesRhoPipeTrks[i] = -1.;
4629 for(Int_t j=0; j<ntracks*2; j++) {
4630 fxOverX0ShieldTrks[j] = -1.;
4631 fxTimesRhoShieldTrks[j] = -1.;
4633 for(Int_t k=0; k<ntracks*6; k++) {
4634 fxOverX0LayerTrks[k] = -1.;
4635 fxTimesRhoLayerTrks[k] = -1.;
4642 //------------------------------------------------------------------------
4643 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4644 //-----------------------------------------------------------------
4645 // Delete LUT for storing material for each prolonged track
4646 //-----------------------------------------------------------------
4647 if(fxOverX0PipeTrks) {
4648 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4650 if(fxOverX0ShieldTrks) {
4651 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4654 if(fxOverX0LayerTrks) {
4655 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4657 if(fxTimesRhoPipeTrks) {
4658 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4660 if(fxTimesRhoShieldTrks) {
4661 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4663 if(fxTimesRhoLayerTrks) {
4664 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4668 //------------------------------------------------------------------------
4669 void AliITStrackerMI::SetForceSkippingOfLayer() {
4670 //-----------------------------------------------------------------
4671 // Check if we are forced to skip layers
4672 // either we set to skip them in RecoParam
4673 // or they were off during data-taking
4674 //-----------------------------------------------------------------
4676 const AliEventInfo *eventInfo = GetEventInfo();
4678 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4679 fForceSkippingOfLayer[l] = 0;
4681 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4685 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4686 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4688 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4689 } else if(l==2 || l==3) {
4690 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4692 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4698 //------------------------------------------------------------------------
4699 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4700 Int_t ilayer,Int_t idet) const {
4701 //-----------------------------------------------------------------
4702 // This method is used to decide whether to allow a prolongation
4703 // without clusters, because we want to skip the layer.
4704 // In this case the return value is > 0:
4705 // return 1: the user requested to skip a layer
4706 // return 2: track outside z acceptance
4707 //-----------------------------------------------------------------
4709 if (ForceSkippingOfLayer(ilayer)) return 1;
4711 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4713 if (idet<0 && // out in z
4714 ilayer>innerLayCanSkip &&
4715 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4716 // check if track will cross SPD outer layer
4717 Double_t phiAtSPD2,zAtSPD2;
4718 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4719 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4721 return 2; // always allow skipping, changed on 05.11.2009
4726 //------------------------------------------------------------------------
4727 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4728 Int_t ilayer,Int_t idet,
4729 Double_t dz,Double_t dy,
4730 Bool_t noClusters) const {
4731 //-----------------------------------------------------------------
4732 // This method is used to decide whether to allow a prolongation
4733 // without clusters, because there is a dead zone in the road.
4734 // In this case the return value is > 0:
4735 // return 1: dead zone at z=0,+-7cm in SPD
4736 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4737 // return 2: all road is "bad" (dead or noisy) from the OCDB
4738 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4739 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4740 //-----------------------------------------------------------------
4742 // check dead zones at z=0,+-7cm in the SPD
4743 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4744 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4745 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4746 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4747 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4748 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4749 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4750 for (Int_t i=0; i<3; i++)
4751 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4752 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4753 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4757 // check bad zones from OCDB
4758 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4760 if (idet<0) return 0;
4762 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4765 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4766 if (ilayer==0 || ilayer==1) { // ---------- SPD
4768 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4770 detSizeFactorX *= 2.;
4771 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4774 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4775 if (detType==2) segm->SetLayer(ilayer+1);
4776 Float_t detSizeX = detSizeFactorX*segm->Dx();
4777 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4779 // check if the road overlaps with bad chips
4781 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4782 Float_t zlocmin = zloc-dz;
4783 Float_t zlocmax = zloc+dz;
4784 Float_t xlocmin = xloc-dy;
4785 Float_t xlocmax = xloc+dy;
4786 Int_t chipsInRoad[100];
4788 // check if road goes out of detector
4789 Bool_t touchNeighbourDet=kFALSE;
4790 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4791 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4792 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4793 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4794 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));
4796 // check if this detector is bad
4798 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4799 if(!touchNeighbourDet) {
4800 return 2; // all detectors in road are bad
4802 return 3; // at least one is bad
4806 if(zlocmin>zlocmax)return 0;
4807 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4808 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4809 if (!nChipsInRoad) return 0;
4811 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4812 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4813 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4814 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4815 if (det.IsChipBad(chipsInRoad[iCh])) {
4823 if(!touchNeighbourDet) {
4824 AliDebug(2,"all bad in road");
4825 return 2; // all chips in road are bad
4827 return 3; // at least a bad chip in road
4832 AliDebug(2,"at least a bad in road");
4833 return 3; // at least a bad chip in road
4837 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4838 || !noClusters) return 0;
4840 // There are no clusters in road: check if there is at least
4841 // a bad SPD pixel or SDD anode or SSD strips on both sides
4843 Int_t idetInITS=idet;
4844 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4846 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4847 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4850 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4854 //------------------------------------------------------------------------
4855 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4856 const AliITStrackMI *track,
4857 Float_t &xloc,Float_t &zloc) const {
4858 //-----------------------------------------------------------------
4859 // Gives position of track in local module ref. frame
4860 //-----------------------------------------------------------------
4865 if(idet<0) return kTRUE; // track out of z acceptance of layer
4867 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4869 Int_t lad = Int_t(idet/ndet) + 1;
4871 Int_t det = idet - (lad-1)*ndet + 1;
4873 Double_t xyzGlob[3],xyzLoc[3];
4875 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4876 // take into account the misalignment: xyz at real detector plane
4877 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4879 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4881 xloc = (Float_t)xyzLoc[0];
4882 zloc = (Float_t)xyzLoc[2];
4886 //------------------------------------------------------------------------
4887 //------------------------------------------------------------------------
4888 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer){
4890 // Method to be optimized further:
4891 // Aim: decide whether a track can be used for PlaneEff evaluation
4892 // the decision is taken based on the track quality at the layer under study
4893 // no information on the clusters on this layer has to be used
4894 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4895 // the cut is done on number of sigmas from the boundaries
4897 // Input: Actual track, layer [0,5] under study
4899 // Return: kTRUE if this is a good track
4901 // it will apply a pre-selection to obtain good quality tracks.
4902 // Here also you will have the possibility to put a control on the
4903 // impact point of the track on the basic block, in order to exclude border regions
4904 // this will be done by calling a proper method of the AliITSPlaneEff class.
4906 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4907 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4909 Int_t index[AliITSgeomTGeo::kNLayers];
4911 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4913 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4914 index[k]=clusters[k];
4918 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4919 AliITSlayer &layer=fgLayers[ilayer];
4920 Double_t r=layer.GetR();
4921 AliITStrackMI tmp(*track);
4923 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4924 Int_t nclout=0; Int_t nclin=0;
4925 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4926 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4927 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4928 // if (tmp.GetClIndex(lay)>=0) nclout++;
4929 if(index[lay]>=0)nclout++;
4931 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4932 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4933 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4934 if (index[lay]>=0) nclin++;
4936 Int_t ncl=nclout+nclin;
4937 Bool_t nextout = kFALSE;
4938 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4939 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4940 Bool_t nextin = kFALSE;
4941 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4942 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4943 // maximum number of missing clusters allowed in outermost layers
4944 if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4946 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4947 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4949 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4950 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4951 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4952 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4956 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4957 Int_t idet=layer.FindDetectorIndex(phi,z);
4958 if(idet<0) { AliInfo(Form("cannot find detector"));
4961 // here check if it has good Chi Square.
4963 //propagate to the intersection with the detector plane
4964 const AliITSdetector &det=layer.GetDetector(idet);
4965 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4969 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4970 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4971 if(key>fPlaneEff->Nblock()) return kFALSE;
4972 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4973 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4975 // DEFINITION OF SEARCH ROAD FOR accepting a track
4977 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4978 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4979 Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXFromBoundaryPlaneEff();
4980 Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZFromBoundaryPlaneEff();
4981 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
4982 // those are precisions in the tracking reference system
4983 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
4984 // Use it also for the module reference system, as it is done for RecPoints
4986 if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigFrmBndPlaneEff()){
4987 if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) dx -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
4990 if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) dz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
4994 // exclude tracks at boundary between detectors
4995 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4996 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4997 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4998 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4999 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5000 if ( (locx-dx < blockXmn+boundaryWidth) ||
5001 (locx+dx > blockXmx-boundaryWidth) ||
5002 (locz-dz < blockZmn+boundaryWidth) ||
5003 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5006 const AliESDEvent *myesd = ((AliESDtrack*)tmp.GetESDtrack())->GetESDEvent();
5008 if (CorrectForPipeMaterial(&tmp,"inward")) {
5009 const AliESDVertex* vtx = myesd->GetVertex();
5010 if(!vtx) return kFALSE;
5011 Double_t ddz[2],cov[3];
5013 if(!tmp.PropagateToDCA(vtx,AliTracker::GetBz(),maxD,ddz,cov)) return kFALSE;
5014 if(TMath::Abs(ddz[0])>=AliITSReconstructor::GetRecoParam()->GetDCACutPlaneEff()) return kFALSE;
5016 Double_t covar[6]; vtx->GetCovMatrix(covar);
5017 Double_t p[2]={tmp.GetParameter()[0]-ddz[0],tmp.GetParameter()[1]-ddz[1]};
5018 Double_t c[3]={covar[2],0.,covar[5]};
5019 Double_t chi2= ((AliESDtrack*)tmp.GetESDtrack())->GetPredictedChi2(p,c);
5020 if (chi2>AliITSReconstructor::GetRecoParam()->GetVertexChi2CutPlaneEff()) return kFALSE; // Use this to cut on chi^2
5022 } else return kFALSE;
5028 //------------------------------------------------------------------------
5029 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
5031 // This Method has to be optimized! For the time-being it uses the same criteria
5032 // as those used in the search of extra clusters for overlapping modules.
5034 // Method Purpose: estabilish whether a track has produced a recpoint or not
5035 // in the layer under study (For Plane efficiency)
5037 // inputs: AliITStrackMI* track (pointer to a usable track)
5039 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5040 // traversed by this very track. In details:
5041 // - if a cluster can be associated to the track then call
5042 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5043 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5046 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5047 AliITSlayer &layer=fgLayers[ilayer];
5048 Double_t r=layer.GetR();
5049 AliITStrackMI tmp(*track);
5053 if (!tmp.GetPhiZat(r,phi,z)) return;
5054 Int_t idet=layer.FindDetectorIndex(phi,z);
5056 if(idet<0) { AliInfo(Form("cannot find detector"));
5060 //propagate to the intersection with the detector plane
5061 const AliITSdetector &det=layer.GetDetector(idet);
5062 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5066 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5068 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5069 TMath::Sqrt(tmp.GetSigmaZ2() +
5070 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5071 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5072 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5073 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5074 TMath::Sqrt(tmp.GetSigmaY2() +
5075 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5076 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5077 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5079 // road in global (rphi,z) [i.e. in tracking ref. system]
5080 Double_t zmin = tmp.GetZ() - dz;
5081 Double_t zmax = tmp.GetZ() + dz;
5082 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5083 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5085 // select clusters in road
5086 layer.SelectClusters(zmin,zmax,ymin,ymax);
5088 // Define criteria for track-cluster association
5089 Double_t msz = tmp.GetSigmaZ2() +
5090 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5091 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5092 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5093 Double_t msy = tmp.GetSigmaY2() +
5094 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5095 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5096 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5097 if (tmp.GetConstrain()) {
5098 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5099 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5101 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5102 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5105 if(AliITSReconstructor::GetRecoParam()->GetSwitchOffStdSearchClusPlaneEff()){
5106 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXSearchClusterPlaneEff();
5107 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZSearchClusterPlaneEff();
5108 Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXSearchClusterPlaneEff();
5109 Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZSearchClusterPlaneEff();
5110 msy = nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
5111 msz = nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
5113 if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigSrhClusPlaneEff()){
5114 if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) msy -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
5117 if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) msz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
5126 if(msz==0 || msy==0){AliWarning("UseTrackForPlaneEff: null search frame"); return;}
5128 msz = 1./msz; // 1/RoadZ^2
5129 msy = 1./msy; // 1/RoadY^2
5131 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5133 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5134 //Double_t tolerance=0.2;
5135 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5136 idetc = cl->GetDetectorIndex();
5137 if(idet!=idetc) continue;
5138 //Int_t ilay = cl->GetLayer();
5140 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5141 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5143 Double_t chi2=tmp.GetPredictedChi2(cl);
5144 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5148 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5150 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5151 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5152 if(key>fPlaneEff->Nblock()) return;
5153 Bool_t found=kFALSE;
5156 while ((cl=layer.GetNextCluster(clidx))!=0) {
5157 idetc = cl->GetDetectorIndex();
5158 if(idet!=idetc) continue;
5159 // here real control to see whether the cluster can be associated to the track.
5160 // cluster not associated to track
5161 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5162 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5163 // calculate track-clusters chi2
5164 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5165 // in particular, the error associated to the cluster
5166 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5168 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5170 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5171 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5172 // track->SetExtraModule(ilayer,idetExtra);
5174 if(!fPlaneEff->UpDatePlaneEff(found,key))
5175 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5177 // this for FO efficiency studies (only for SPD) //
5178 UInt_t keyFO=999999;
5179 Bool_t foundFO=kFALSE;
5180 if(ilayer<2){ //ONLY SPD layers for FastOr studies
5181 TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5182 Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5183 if(!fSPDChipIntPlaneEff[key]){
5184 AliITSPlaneEffSPD spd;
5185 keyFO = spd.SwitchChipKeyNumbering(key);
5186 if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5187 keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5188 if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5189 AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5192 if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5193 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5199 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5200 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5201 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5202 Int_t cltype[2]={-999,-999};
5205 Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5209 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5210 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5213 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5214 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5215 cltype[0]=layer.GetCluster(ci)->GetNy();
5216 cltype[1]=layer.GetCluster(ci)->GetNz();
5218 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5219 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5220 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5221 // It is computed properly by calling the method
5222 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5224 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5225 //if (tmp.PropagateTo(x,0.,0.)) {
5226 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5227 AliCluster c(*layer.GetCluster(ci));
5228 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5229 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5230 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5231 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5232 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5235 // Compute the angles between the track and the module
5236 // compute the angle "in phi direction", i.e. the angle in the transverse plane
5237 // between the normal to the module and the projection (in the transverse plane) of the
5239 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5240 Float_t tgl = tmp.GetTgl();
5241 Float_t phitr = tmp.GetSnp();
5242 phitr = TMath::ASin(phitr);
5243 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5245 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5246 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5248 alpha = tmp.GetAlpha();
5249 Double_t phiglob = alpha+phitr;
5251 p[0] = TMath::Cos(phiglob);
5252 p[1] = TMath::Sin(phiglob);
5254 TVector3 pvec(p[0],p[1],p[2]);
5255 TVector3 normvec(rot[1],rot[4],rot[7]);
5256 Double_t angle = pvec.Angle(normvec);
5258 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5259 angle *= 180./TMath::Pi();
5262 TVector3 pt(p[0],p[1],0);
5263 TVector3 normt(rot[1],rot[4],0);
5264 Double_t anglet = pt.Angle(normt);
5266 Double_t phiPt = TMath::ATan2(p[1],p[0]);
5267 if(phiPt<0)phiPt+=2.*TMath::Pi();
5268 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5269 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5270 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5271 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
5272 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5273 anglet *= 180./TMath::Pi();
5275 angleModTrack[2]=(Float_t) angle;
5276 angleModTrack[0]=(Float_t) anglet;
5277 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5278 angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5279 angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5280 angleModTrack[1]*=180./TMath::Pi(); // in degree
5282 fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
5284 // For FO efficiency studies of SPD
5285 if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
5287 if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5291 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5293 // find the MC cluster for the label
5294 return fgLayers[lr].FindClusterForLabel(label,store);
5298 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5300 // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5301 strncpy(patt,"......",6);
5303 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5306 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5307 Int_t cindex = track->GetClusterIndex(i);
5308 Int_t l=(cindex & 0xf0000000) >> 28;
5309 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5311 for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5312 patt[l] = isWrong ? 'f':'c';
5318 //------------------------------------------------------------------------
5319 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5321 //--------------------------------------------------------------------
5324 for (int ic=0;ic<fN;ic++) {
5325 const AliITSRecPoint *cl = GetCluster(ic);
5326 for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5328 if (store) store[nfound] = ic;