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>
38 #include "AliGeomManager.h"
39 #include "AliITSPlaneEff.h"
40 #include "AliTrackPointArray.h"
41 #include "AliESDEvent.h"
42 #include "AliESDV0Params.h"
43 #include "AliESDtrack.h"
45 #include "AliITSChannelStatus.h"
46 #include "AliITSDetTypeRec.h"
47 #include "AliITSRecPoint.h"
48 #include "AliITSRecPointContainer.h"
49 #include "AliITSgeomTGeo.h"
50 #include "AliITSReconstructor.h"
51 #include "AliITSClusterParam.h"
52 #include "AliITSsegmentation.h"
53 #include "AliITSCalibration.h"
54 #include "AliITSPlaneEffSPD.h"
55 #include "AliITSPlaneEffSDD.h"
56 #include "AliITSPlaneEffSSD.h"
57 #include "AliITSV0Finder.h"
58 #include "AliITStrackerMI.h"
59 #include "AliMathBase.h"
62 ClassImp(AliITStrackerMI)
64 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
66 AliITStrackerMI::AliITStrackerMI():AliTracker(),
76 fLastLayerToTrackTo(0),
79 fTrackingPhase("Default"),
83 fSelectBestMIP03(kFALSE),
84 fUseImproveKalman(kFALSE),
88 fxTimesRhoPipeTrks(0),
89 fxOverX0ShieldTrks(0),
90 fxTimesRhoShieldTrks(0),
92 fxTimesRhoLayerTrks(0),
100 //Default constructor
102 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
104 fxOverX0Shield[i]=-1.;
105 fxTimesRhoShield[i]=-1.;
108 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
109 fOriginal.SetOwner();
110 for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
111 for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
112 fITSPid=new AliITSPIDResponse();
115 //------------------------------------------------------------------------
116 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
117 fI(AliITSgeomTGeo::GetNLayers()),
126 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
129 fTrackingPhase("Default"),
133 fSelectBestMIP03(kFALSE),
134 fUseImproveKalman(kFALSE),
138 fxTimesRhoPipeTrks(0),
139 fxOverX0ShieldTrks(0),
140 fxTimesRhoShieldTrks(0),
141 fxOverX0LayerTrks(0),
142 fxTimesRhoLayerTrks(0),
144 fITSChannelStatus(0),
148 //--------------------------------------------------------------------
149 //This is the AliITStrackerMI constructor
150 //--------------------------------------------------------------------
152 AliWarning("\"geom\" is actually a dummy argument !");
155 fOriginal.SetOwner();
159 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
160 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
161 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
163 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
164 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
165 Double_t poff=TMath::ATan2(y,x);
168 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
169 Double_t r=TMath::Sqrt(x*x + y*y);
171 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
172 r += TMath::Sqrt(x*x + y*y);
173 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
174 r += TMath::Sqrt(x*x + y*y);
175 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
176 r += TMath::Sqrt(x*x + y*y);
179 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
181 for (Int_t j=1; j<nlad+1; j++) {
182 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
183 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
184 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
186 Double_t txyz[3]={0.};
187 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
188 m.LocalToMaster(txyz,xyz);
189 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
190 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
192 if (phi<0) phi+=TMath::TwoPi();
193 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
195 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
196 new(&det) AliITSdetector(r,phi);
197 // compute the real radius (with misalignment)
198 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
200 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
201 mmisal.LocalToMaster(txyz,xyz);
202 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
203 det.SetRmisal(rmisal);
205 } // end loop on detectors
206 } // end loop on ladders
207 fForceSkippingOfLayer[i-1] = 0;
208 } // end loop on layers
211 fI=AliITSgeomTGeo::GetNLayers();
214 fConstraint[0]=1; fConstraint[1]=0;
216 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
217 AliITSReconstructor::GetRecoParam()->GetYVdef(),
218 AliITSReconstructor::GetRecoParam()->GetZVdef()};
219 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
220 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
221 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
222 SetVertex(xyzVtx,ersVtx);
224 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
225 for (Int_t i=0;i<100000;i++){
226 fBestTrackIndex[i]=0;
229 // store positions of centre of SPD modules (in z)
230 // The convetion is that fSPDdetzcentre is ordered from -z to +z
232 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
233 if (tr[2]<0) { // old geom (up to v5asymmPPR)
234 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
235 fSPDdetzcentre[0] = tr[2];
236 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
237 fSPDdetzcentre[1] = tr[2];
238 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
239 fSPDdetzcentre[2] = tr[2];
240 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
241 fSPDdetzcentre[3] = tr[2];
242 } else { // new geom (from v11Hybrid)
243 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
244 fSPDdetzcentre[0] = tr[2];
245 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
246 fSPDdetzcentre[1] = tr[2];
247 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
248 fSPDdetzcentre[2] = tr[2];
249 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
250 fSPDdetzcentre[3] = tr[2];
253 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
254 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
255 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
259 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
260 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
262 if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
263 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
265 // only for plane efficiency evaluation
266 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
267 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
268 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
269 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
270 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
271 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
272 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
273 else fPlaneEff = new AliITSPlaneEffSSD();
274 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
275 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
276 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
280 fSelectBestMIP03 = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
281 fFlagFakes = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
282 fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
284 fITSPid=new AliITSPIDResponse();
287 //------------------------------------------------------------------------
288 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
290 fBestTrack(tracker.fBestTrack),
291 fTrackToFollow(tracker.fTrackToFollow),
292 fTrackHypothesys(tracker.fTrackHypothesys),
293 fBestHypothesys(tracker.fBestHypothesys),
294 fOriginal(tracker.fOriginal),
295 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
296 fPass(tracker.fPass),
297 fAfterV0(tracker.fAfterV0),
298 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
299 fCoefficients(tracker.fCoefficients),
301 fTrackingPhase(tracker.fTrackingPhase),
302 fUseTGeo(tracker.fUseTGeo),
303 fNtracks(tracker.fNtracks),
304 fFlagFakes(tracker.fFlagFakes),
305 fSelectBestMIP03(tracker.fSelectBestMIP03),
306 fxOverX0Pipe(tracker.fxOverX0Pipe),
307 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
309 fxTimesRhoPipeTrks(0),
310 fxOverX0ShieldTrks(0),
311 fxTimesRhoShieldTrks(0),
312 fxOverX0LayerTrks(0),
313 fxTimesRhoLayerTrks(0),
314 fDebugStreamer(tracker.fDebugStreamer),
315 fITSChannelStatus(tracker.fITSChannelStatus),
316 fkDetTypeRec(tracker.fkDetTypeRec),
317 fPlaneEff(tracker.fPlaneEff) {
319 fOriginal.SetOwner();
322 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
325 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
326 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
329 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
330 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
334 //------------------------------------------------------------------------
335 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
336 //Assignment operator
337 this->~AliITStrackerMI();
338 new(this) AliITStrackerMI(tracker);
342 //------------------------------------------------------------------------
343 AliITStrackerMI::~AliITStrackerMI()
348 if (fCoefficients) delete [] fCoefficients;
349 DeleteTrksMaterialLUT();
350 if (fDebugStreamer) {
351 //fDebugStreamer->Close();
352 delete fDebugStreamer;
354 if(fITSChannelStatus) delete fITSChannelStatus;
355 if(fPlaneEff) delete fPlaneEff;
356 if(fITSPid) delete fITSPid;
359 //------------------------------------------------------------------------
360 void AliITStrackerMI::ReadBadFromDetTypeRec() {
361 //--------------------------------------------------------------------
362 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
364 //--------------------------------------------------------------------
366 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
368 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
370 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
373 if(fITSChannelStatus) delete fITSChannelStatus;
374 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
376 // ITS detectors and chips
377 Int_t i=0,j=0,k=0,ndet=0;
378 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
379 Int_t nBadDetsPerLayer=0;
380 ndet=AliITSgeomTGeo::GetNDetectors(i);
381 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
382 for (k=1; k<ndet+1; k++) {
383 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
384 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
385 if(det.IsBad()) {nBadDetsPerLayer++;}
386 } // end loop on detectors
387 } // end loop on ladders
388 AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
389 } // end loop on layers
393 //------------------------------------------------------------------------
394 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
395 //--------------------------------------------------------------------
396 //This function loads ITS clusters
397 //--------------------------------------------------------------------
399 TClonesArray *clusters = NULL;
400 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
401 clusters=rpcont->FetchClusters(0,cTree);
402 if(!clusters) return 1;
404 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
405 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
408 Int_t i=0,j=0,ndet=0;
410 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
411 ndet=fgLayers[i].GetNdetectors();
412 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
413 for (; j<jmax; j++) {
414 // if (!cTree->GetEvent(j)) continue;
415 clusters = rpcont->UncheckedGetClusters(j);
416 if(!clusters)continue;
417 Int_t ncl=clusters->GetEntriesFast();
418 SignDeltas(clusters,GetZ());
421 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
422 detector=c->GetDetectorIndex();
424 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
426 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
428 AliWarning(Form("Too many clusters on layer %d!",i));
433 // add dead zone "virtual" cluster in SPD, if there is a cluster within
434 // zwindow cm from the dead zone
435 // This method assumes that fSPDdetzcentre is ordered from -z to +z
436 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
437 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
438 Int_t lab[4] = {0,0,0,detector};
439 Int_t info[3] = {0,0,i};
440 Float_t q = 0.; // this identifies virtual clusters
441 Float_t hit[6] = {xdead,
443 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
444 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
447 Bool_t local = kTRUE;
448 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
449 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
450 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
451 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
452 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
453 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
454 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
455 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
456 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
457 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
458 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
459 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
460 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
461 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
462 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
463 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
464 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
465 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
466 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
468 } // "virtual" clusters in SPD
472 fgLayers[i].ResetRoad(); //road defined by the cluster density
473 fgLayers[i].SortClusters();
476 // check whether we have to skip some layers
477 SetForceSkippingOfLayer();
481 //------------------------------------------------------------------------
482 void AliITStrackerMI::UnloadClusters() {
483 //--------------------------------------------------------------------
484 //This function unloads ITS clusters
485 //--------------------------------------------------------------------
486 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
488 //------------------------------------------------------------------------
489 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
490 //--------------------------------------------------------------------
491 // Publishes all pointers to clusters known to the tracker into the
492 // passed object array.
493 // The ownership is not transfered - the caller is not expected to delete
495 //--------------------------------------------------------------------
497 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
498 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
499 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
506 //------------------------------------------------------------------------
507 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
508 //--------------------------------------------------------------------
509 // Correction for the material between the TPC and the ITS
510 //--------------------------------------------------------------------
511 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
512 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
513 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
514 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
515 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
516 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
517 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
518 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
520 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
526 //------------------------------------------------------------------------
527 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
528 //--------------------------------------------------------------------
529 // This functions reconstructs ITS tracks
530 // The clusters must be already loaded !
531 //--------------------------------------------------------------------
533 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
535 fTrackingPhase="Clusters2Tracks";
538 fSelectBestMIP03 = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
539 fFlagFakes = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
540 fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
542 TObjArray itsTracks(15000);
544 fEsd = event; // store pointer to the esd
546 // temporary (for cosmics)
547 if(event->GetVertex()) {
548 TString title = event->GetVertex()->GetTitle();
549 if(title.Contains("cosmics")) {
550 Double_t xyz[3]={GetX(),GetY(),GetZ()};
551 Double_t exyz[3]={0.1,0.1,0.1};
557 {/* Read ESD tracks */
558 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
559 Int_t nentr=event->GetNumberOfTracks();
561 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
563 AliESDtrack *esd=event->GetTrack(nentr);
564 // ---- for debugging:
565 //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); }
567 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
568 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
569 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
570 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
571 AliITStrackMI *t = new AliITStrackMI(*esd);
572 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
573 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
576 // look at the ESD mass hypothesys !
577 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
579 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
581 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
582 //track - can be V0 according to TPC
584 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
588 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
592 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
596 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
601 t->SetReconstructed(kFALSE);
602 itsTracks.AddLast(t);
603 fOriginal.AddLast(t);
605 } /* End Read ESD tracks */
609 Int_t nentr=itsTracks.GetEntriesFast();
610 fTrackHypothesys.Expand(nentr);
611 fBestHypothesys.Expand(nentr);
612 MakeCoefficients(nentr);
613 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
615 // THE TWO TRACKING PASSES
616 for (fPass=0; fPass<2; fPass++) {
617 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
618 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
619 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
620 if (t==0) continue; //this track has been already tracked
621 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
622 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
623 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
624 if (fConstraint[fPass]) {
625 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
626 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
629 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
630 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
632 ResetTrackToFollow(*t);
635 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
638 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
640 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
641 if (!besttrack) continue;
642 besttrack->SetLabel(tpcLabel);
643 // besttrack->CookdEdx();
645 besttrack->SetFakeRatio(1.);
646 CookLabel(besttrack,0.); //For comparison only
647 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
648 t->SetWinner(besttrack);
650 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
652 t->SetReconstructed(kTRUE);
654 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
656 GetBestHypothesysMIP(itsTracks);
657 } // end loop on the two tracking passes
659 if (fFlagFakes) FlagFakes(itsTracks);
661 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
662 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
667 Int_t entries = fTrackHypothesys.GetEntriesFast();
668 for (Int_t ientry=0; ientry<entries; ientry++) {
669 TObjArray * array =(TObjArray*)fTrackHypothesys.At(ientry);
670 if (array) array->Delete();
671 delete fTrackHypothesys.RemoveAt(ientry);
674 fTrackHypothesys.Delete();
675 entries = fBestHypothesys.GetEntriesFast();
676 for (Int_t ientry=0; ientry<entries; ientry++) {
677 TObjArray * array =(TObjArray*)fBestHypothesys.At(ientry);
678 if (array) array->Delete();
679 delete fBestHypothesys.RemoveAt(ientry);
681 fBestHypothesys.Delete();
683 delete [] fCoefficients;
685 DeleteTrksMaterialLUT();
687 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
689 fTrackingPhase="Default";
693 //------------------------------------------------------------------------
694 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
695 //--------------------------------------------------------------------
696 // This functions propagates reconstructed ITS tracks back
697 // The clusters must be loaded !
698 //--------------------------------------------------------------------
699 fTrackingPhase="PropagateBack";
700 Int_t nentr=event->GetNumberOfTracks();
701 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
704 for (Int_t i=0; i<nentr; i++) {
705 AliESDtrack *esd=event->GetTrack(i);
707 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
708 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
710 AliITStrackMI *t = new AliITStrackMI(*esd);
712 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
714 ResetTrackToFollow(*t);
717 // propagate to vertex [SR, GSI 17.02.2003]
718 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
719 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
720 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
721 fTrackToFollow.StartTimeIntegral();
722 // from vertex to outside pipe
723 CorrectForPipeMaterial(&fTrackToFollow,"outward");
725 // Start time integral and add distance from current position to vertex
726 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
727 fTrackToFollow.GetXYZ(xyzTrk);
729 for (Int_t icoord=0; icoord<3; icoord++) {
730 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
733 fTrackToFollow.StartTimeIntegral();
734 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
736 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
737 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
738 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
742 fTrackToFollow.SetLabel(t->GetLabel());
743 //fTrackToFollow.CookdEdx();
744 CookLabel(&fTrackToFollow,0.); //For comparison only
745 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
746 //UseClusters(&fTrackToFollow);
752 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
754 fTrackingPhase="Default";
758 //------------------------------------------------------------------------
759 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
760 //--------------------------------------------------------------------
761 // This functions refits ITS tracks using the
762 // "inward propagated" TPC tracks
763 // The clusters must be loaded !
764 //--------------------------------------------------------------------
765 fTrackingPhase="RefitInward";
767 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
769 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
770 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
772 Int_t nentr=event->GetNumberOfTracks();
773 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
776 for (Int_t i=0; i<nentr; i++) {
777 AliESDtrack *esd=event->GetTrack(i);
779 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
780 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
781 if (esd->GetStatus()&AliESDtrack::kTPCout)
782 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
784 AliITStrackMI *t = new AliITStrackMI(*esd);
786 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
787 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
792 ResetTrackToFollow(*t);
793 fTrackToFollow.ResetClusters();
795 // ITS standalone tracks
796 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
797 fTrackToFollow.ResetCovariance(10.);
798 // protection for loopers that can have parameters screwed up
799 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
800 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
807 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
808 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
810 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
811 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
812 AliDebug(2," refit OK");
813 fTrackToFollow.SetLabel(t->GetLabel());
814 // fTrackToFollow.CookdEdx();
815 CookdEdx(&fTrackToFollow);
817 CookLabel(&fTrackToFollow,0.0); //For comparison only
820 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
821 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
822 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
823 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
824 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
825 Double_t r[3]={0.,0.,0.};
827 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
834 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
836 fTrackingPhase="Default";
840 //------------------------------------------------------------------------
841 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
842 //--------------------------------------------------------------------
843 // Return pointer to a given cluster
844 //--------------------------------------------------------------------
845 Int_t l=(index & 0xf0000000) >> 28;
846 Int_t c=(index & 0x0fffffff) >> 00;
847 return fgLayers[l].GetCluster(c);
849 //------------------------------------------------------------------------
850 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
851 //--------------------------------------------------------------------
852 // Get track space point with index i
853 //--------------------------------------------------------------------
855 Int_t l=(index & 0xf0000000) >> 28;
856 Int_t c=(index & 0x0fffffff) >> 00;
857 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
858 Int_t idet = cl->GetDetectorIndex();
862 cl->GetGlobalXYZ(xyz);
863 cl->GetGlobalCov(cov);
865 p.SetCharge(cl->GetQ());
866 p.SetDriftTime(cl->GetDriftTime());
867 p.SetChargeRatio(cl->GetChargeRatio());
868 p.SetClusterType(cl->GetClusterType());
869 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
872 iLayer = AliGeomManager::kSPD1;
875 iLayer = AliGeomManager::kSPD2;
878 iLayer = AliGeomManager::kSDD1;
881 iLayer = AliGeomManager::kSDD2;
884 iLayer = AliGeomManager::kSSD1;
887 iLayer = AliGeomManager::kSSD2;
890 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
893 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
894 p.SetVolumeID((UShort_t)volid);
897 //------------------------------------------------------------------------
898 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
899 AliTrackPoint& p, const AliESDtrack *t) {
900 //--------------------------------------------------------------------
901 // Get track space point with index i
902 // (assign error estimated during the tracking)
903 //--------------------------------------------------------------------
905 Int_t l=(index & 0xf0000000) >> 28;
906 Int_t c=(index & 0x0fffffff) >> 00;
907 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
908 Int_t idet = cl->GetDetectorIndex();
910 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
912 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
914 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
915 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
916 Double_t alpha = t->GetAlpha();
917 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
918 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
919 phi += alpha-det.GetPhi();
920 Float_t tgphi = TMath::Tan(phi);
922 Float_t tgl = t->GetTgl(); // tgl about const along track
923 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
925 Float_t errtrky,errtrkz,covyz;
926 Bool_t addMisalErr=kFALSE;
927 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
931 cl->GetGlobalXYZ(xyz);
932 // cl->GetGlobalCov(cov);
933 Float_t pos[3] = {0.,0.,0.};
934 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
935 tmpcl.GetGlobalCov(cov);
938 p.SetCharge(cl->GetQ());
939 p.SetDriftTime(cl->GetDriftTime());
940 p.SetChargeRatio(cl->GetChargeRatio());
941 p.SetClusterType(cl->GetClusterType());
943 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
946 iLayer = AliGeomManager::kSPD1;
949 iLayer = AliGeomManager::kSPD2;
952 iLayer = AliGeomManager::kSDD1;
955 iLayer = AliGeomManager::kSDD2;
958 iLayer = AliGeomManager::kSSD1;
961 iLayer = AliGeomManager::kSSD2;
964 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
967 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
969 p.SetVolumeID((UShort_t)volid);
972 //------------------------------------------------------------------------
973 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
975 //--------------------------------------------------------------------
976 // Follow prolongation tree
977 //--------------------------------------------------------------------
979 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
980 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
983 AliESDtrack * esd = otrack->GetESDtrack();
984 if (esd->GetV0Index(0)>0) {
985 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
986 // mapping of ESD track is different as ITS track in Containers
987 // Need something more stable
988 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
989 for (Int_t i=0;i<3;i++){
990 Int_t index = esd->GetV0Index(i);
992 AliESDv0 * vertex = fEsd->GetV0(index);
993 if (vertex->GetStatus()<0) continue; // rejected V0
995 if (esd->GetSign()>0) {
996 vertex->SetIndex(0,esdindex);
998 vertex->SetIndex(1,esdindex);
1002 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1004 bestarray = new TObjArray(5);
1005 bestarray->SetOwner();
1006 fBestHypothesys.AddAt(bestarray,esdindex);
1010 //setup tree of the prolongations
1012 const int kMaxTr = 100; //RS
1013 static AliITStrackMI tracks[7][kMaxTr];
1014 AliITStrackMI *currenttrack;
1015 static AliITStrackMI currenttrack1;
1016 static AliITStrackMI currenttrack2;
1017 static AliITStrackMI backuptrack;
1019 Int_t nindexes[7][kMaxTr];
1020 Float_t normalizedchi2[kMaxTr];
1021 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1022 otrack->SetNSkipped(0);
1023 new (&(tracks[6][0])) AliITStrackMI(*otrack);
1025 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1026 Int_t modstatus = 1; // found
1030 // follow prolongations
1031 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1032 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1035 AliITSlayer &layer=fgLayers[ilayer];
1036 Double_t r = layer.GetR();
1042 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1043 // printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1045 if (ntracks[ilayer]>=kMaxTr) break;
1046 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1047 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1048 if (ntracks[ilayer]>15+ilayer){
1049 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1050 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1053 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1055 // material between SSD and SDD, SDD and SPD
1057 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1059 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1063 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1064 Int_t idet=layer.FindDetectorIndex(phi,z);
1066 Double_t trackGlobXYZ1[3];
1067 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1069 // Get the budget to the primary vertex for the current track being prolonged
1070 Double_t budgetToPrimVertex = 0;
1071 double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1074 if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1075 else budgetToPrimVertex = GetEffectiveThickness();
1077 // check if we allow a prolongation without point
1078 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1080 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1081 // propagate to the layer radius
1082 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1083 if(!vtrack->Propagate(xToGo)) continue;
1084 // apply correction for material of the current layer
1085 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1086 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1087 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1088 vtrack->SetClIndex(ilayer,-1);
1089 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1090 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1091 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1093 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1095 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1096 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1102 // track outside layer acceptance in z
1103 if (idet<0) continue;
1105 //propagate to the intersection with the detector plane
1106 const AliITSdetector &det=layer.GetDetector(idet);
1107 new(¤ttrack2) AliITStrackMI(currenttrack1);
1108 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1109 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1110 currenttrack1.SetDetectorIndex(idet);
1111 currenttrack2.SetDetectorIndex(idet);
1112 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1115 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1117 // road in global (rphi,z) [i.e. in tracking ref. system]
1118 Double_t zmin,zmax,ymin,ymax;
1119 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1121 // select clusters in road
1122 layer.SelectClusters(zmin,zmax,ymin,ymax);
1123 //********************
1125 // Define criteria for track-cluster association
1126 Double_t msz = currenttrack1.GetSigmaZ2() +
1127 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1128 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1129 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1130 Double_t msy = currenttrack1.GetSigmaY2() +
1131 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1132 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1133 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1135 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1136 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1138 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1139 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1141 msz = 1./msz; // 1/RoadZ^2
1142 msy = 1./msy; // 1/RoadY^2
1146 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1148 const AliITSRecPoint *cl=0;
1150 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1151 Bool_t deadzoneSPD=kFALSE;
1152 currenttrack = ¤ttrack1;
1154 // check if the road contains a dead zone
1155 Bool_t noClusters = kFALSE;
1156 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1157 if (noClusters) AliDebug(2,"no clusters in road");
1158 Double_t dz=0.5*(zmax-zmin);
1159 Double_t dy=0.5*(ymax-ymin);
1160 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1161 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1162 // create a prolongation without clusters (check also if there are no clusters in the road)
1165 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1166 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1167 updatetrack->SetClIndex(ilayer,-1);
1169 modstatus = 5; // no cls in road
1170 } else if (dead==1) {
1171 modstatus = 7; // holes in z in SPD
1172 } else if (dead==2 || dead==3 || dead==4) {
1173 modstatus = 2; // dead from OCDB
1175 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1176 // apply correction for material of the current layer
1177 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1178 if (constrain) { // apply vertex constrain
1179 updatetrack->SetConstrain(constrain);
1180 Bool_t isPrim = kTRUE;
1181 if (ilayer<4) { // check that it's close to the vertex
1182 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1183 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1184 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1185 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1186 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1188 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1190 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1191 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1194 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1196 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1197 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1199 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1200 updatetrack->SetDeadZoneProbability(ilayer,1.);
1201 } else if (dead==4) { // at least a single dead channel from OCDB
1202 updatetrack->SetDeadZoneProbability(ilayer,0.);
1209 // loop over clusters in the road
1210 while ((cl=layer.GetNextCluster(clidx))!=0) {
1211 if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters
1212 Bool_t changedet =kFALSE;
1213 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1214 Int_t idetc=cl->GetDetectorIndex();
1216 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1217 // take into account misalignment (bring track to real detector plane)
1218 Double_t xTrOrig = currenttrack->GetX();
1219 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1220 // a first cut on track-cluster distance
1221 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1222 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1223 { // cluster not associated to track
1224 AliDebug(2,"not associated");
1225 // MvL: added here as well
1226 // bring track back to ideal detector plane
1227 currenttrack->Propagate(xTrOrig);
1230 // bring track back to ideal detector plane
1231 if (!currenttrack->Propagate(xTrOrig)) continue;
1232 } else { // have to move track to cluster's detector
1233 const AliITSdetector &detc=layer.GetDetector(idetc);
1234 // a first cut on track-cluster distance
1236 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1237 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1238 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1239 continue; // cluster not associated to track
1241 new (&backuptrack) AliITStrackMI(currenttrack2);
1243 currenttrack =¤ttrack2;
1244 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1245 new (currenttrack) AliITStrackMI(backuptrack);
1249 currenttrack->SetDetectorIndex(idetc);
1250 // Get again the budget to the primary vertex
1251 // for the current track being prolonged, if had to change detector
1252 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1255 // calculate track-clusters chi2
1256 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1258 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1259 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1260 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1261 if (ntracks[ilayer]>=kMaxTr) continue;
1262 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1263 updatetrack->SetClIndex(ilayer,-1);
1264 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1266 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1267 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1268 AliDebug(2,"update failed");
1271 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1272 modstatus = 1; // found
1273 } else { // virtual cluster in dead zone
1274 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1275 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1276 modstatus = 7; // holes in z in SPD
1280 Float_t xlocnewdet,zlocnewdet;
1281 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1282 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1285 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1287 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1289 // apply correction for material of the current layer
1290 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1292 if (constrain) { // apply vertex constrain
1293 updatetrack->SetConstrain(constrain);
1294 Bool_t isPrim = kTRUE;
1295 if (ilayer<4) { // check that it's close to the vertex
1296 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1297 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1298 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1299 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1300 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1302 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1304 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1305 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1307 } //apply vertex constrain
1309 } // create new hypothesis
1311 AliDebug(2,"chi2 too large");
1314 } // loop over possible prolongations
1316 // allow one prolongation without clusters
1317 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1318 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1319 // apply correction for material of the current layer
1320 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1321 vtrack->SetClIndex(ilayer,-1);
1322 modstatus = 3; // skipped
1323 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1324 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1326 vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1327 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1329 vtrack->IncrementNSkipped();
1334 } // loop over tracks in layer ilayer+1
1336 //loop over track candidates for the current layer
1342 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1343 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1344 if (normalizedchi2[itrack] <
1345 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1349 if (constrain) { // constrain
1350 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1352 } else { // !constrain
1353 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1358 // sort tracks by increasing normalized chi2
1359 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1360 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1361 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1362 // if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1363 if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9);
1364 } // end loop over layers
1368 // Now select tracks to be kept
1370 Int_t max = constrain ? 20 : 5;
1372 // tracks that reach layer 0 (SPD inner)
1373 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1374 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1375 if (track.GetNumberOfClusters()<2) continue;
1376 if (!constrain && track.GetNormChi2(0) >
1377 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1380 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1383 // tracks that reach layer 1 (SPD outer)
1384 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1385 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1386 if (track.GetNumberOfClusters()<4) continue;
1387 if (!constrain && track.GetNormChi2(1) >
1388 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1389 if (constrain) track.IncrementNSkipped();
1391 track.SetD(0,track.GetD(GetX(),GetY()));
1392 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1393 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1394 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1397 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1400 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1402 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1403 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1404 if (track.GetNumberOfClusters()<3) continue;
1405 if (track.GetNormChi2(2) >
1406 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1407 track.SetD(0,track.GetD(GetX(),GetY()));
1408 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1409 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1410 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1412 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1418 // register best track of each layer - important for V0 finder
1420 for (Int_t ilayer=0;ilayer<5;ilayer++){
1421 if (ntracks[ilayer]==0) continue;
1422 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1423 if (track.GetNumberOfClusters()<1) continue;
1424 CookLabel(&track,0);
1425 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1429 // update TPC V0 information
1431 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1432 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1433 for (Int_t i=0;i<3;i++){
1434 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1435 if (index==0) break;
1436 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1437 if (vertex->GetStatus()<0) continue; // rejected V0
1439 if (otrack->GetSign()>0) {
1440 vertex->SetIndex(0,esdindex);
1443 vertex->SetIndex(1,esdindex);
1445 //find nearest layer with track info
1446 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1447 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1448 Int_t nearest = nearestold;
1449 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1450 if (ntracks[nearest]==0){
1455 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1456 if (nearestold<5&&nearest<5){
1457 Bool_t accept = track.GetNormChi2(nearest)<10;
1459 if (track.GetSign()>0) {
1460 vertex->SetParamP(track);
1461 vertex->Update(fprimvertex);
1462 //vertex->SetIndex(0,track.fESDtrack->GetID());
1463 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1465 vertex->SetParamN(track);
1466 vertex->Update(fprimvertex);
1467 //vertex->SetIndex(1,track.fESDtrack->GetID());
1468 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1470 vertex->SetStatus(vertex->GetStatus()+1);
1472 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1479 //------------------------------------------------------------------------
1480 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1482 //--------------------------------------------------------------------
1485 return fgLayers[layer];
1487 //------------------------------------------------------------------------
1488 AliITStrackerMI::AliITSlayer::AliITSlayer():
1518 //--------------------------------------------------------------------
1519 //default AliITSlayer constructor
1520 //--------------------------------------------------------------------
1521 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1522 fClusterWeight[i]=0;
1523 fClusterTracks[0][i]=-1;
1524 fClusterTracks[1][i]=-1;
1525 fClusterTracks[2][i]=-1;
1526 fClusterTracks[3][i]=-1;
1533 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1534 for (Int_t j1=0; j1<6; j1++) {
1535 fClusters5[j1][j]=0;
1536 fClusterIndex5[j1][j]=-1;
1545 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1546 for (Int_t j1=0; j1<11; j1++) {
1547 fClusters10[j1][j]=0;
1548 fClusterIndex10[j1][j]=-1;
1557 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1558 for (Int_t j1=0; j1<21; j1++) {
1559 fClusters20[j1][j]=0;
1560 fClusterIndex20[j1][j]=-1;
1568 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1573 //------------------------------------------------------------------------
1574 AliITStrackerMI::AliITSlayer::
1575 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1604 //--------------------------------------------------------------------
1605 //main AliITSlayer constructor
1606 //--------------------------------------------------------------------
1607 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1608 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1610 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1611 fClusterWeight[i]=0;
1612 fClusterTracks[0][i]=-1;
1613 fClusterTracks[1][i]=-1;
1614 fClusterTracks[2][i]=-1;
1615 fClusterTracks[3][i]=-1;
1623 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1624 for (Int_t j1=0; j1<6; j1++) {
1625 fClusters5[j1][j]=0;
1626 fClusterIndex5[j1][j]=-1;
1635 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1636 for (Int_t j1=0; j1<11; j1++) {
1637 fClusters10[j1][j]=0;
1638 fClusterIndex10[j1][j]=-1;
1647 for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1648 for (Int_t j1=0; j1<21; j1++) {
1649 fClusters20[j1][j]=0;
1650 fClusterIndex20[j1][j]=-1;
1658 for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1664 //------------------------------------------------------------------------
1665 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1667 fPhiOffset(layer.fPhiOffset),
1668 fNladders(layer.fNladders),
1669 fZOffset(layer.fZOffset),
1670 fNdetectors(layer.fNdetectors),
1671 fDetectors(layer.fDetectors),
1676 fClustersCs(layer.fClustersCs),
1677 fClusterIndexCs(layer.fClusterIndexCs),
1681 fCurrentSlice(layer.fCurrentSlice),
1689 fAccepted(layer.fAccepted),
1691 fMaxSigmaClY(layer.fMaxSigmaClY),
1692 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1693 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1698 //------------------------------------------------------------------------
1699 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1700 //--------------------------------------------------------------------
1701 // AliITSlayer destructor
1702 //--------------------------------------------------------------------
1703 delete [] fDetectors;
1704 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1705 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1706 fClusterWeight[i]=0;
1707 fClusterTracks[0][i]=-1;
1708 fClusterTracks[1][i]=-1;
1709 fClusterTracks[2][i]=-1;
1710 fClusterTracks[3][i]=-1;
1713 //------------------------------------------------------------------------
1714 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1715 //--------------------------------------------------------------------
1716 // This function removes loaded clusters
1717 //--------------------------------------------------------------------
1718 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1719 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1720 fClusterWeight[i]=0;
1721 fClusterTracks[0][i]=-1;
1722 fClusterTracks[1][i]=-1;
1723 fClusterTracks[2][i]=-1;
1724 fClusterTracks[3][i]=-1;
1730 //------------------------------------------------------------------------
1731 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1732 //--------------------------------------------------------------------
1733 // This function reset weights of the clusters
1734 //--------------------------------------------------------------------
1735 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1736 fClusterWeight[i]=0;
1737 fClusterTracks[0][i]=-1;
1738 fClusterTracks[1][i]=-1;
1739 fClusterTracks[2][i]=-1;
1740 fClusterTracks[3][i]=-1;
1742 for (Int_t i=0; i<fN;i++) {
1743 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1744 if (cl&&cl->IsUsed()) cl->Use();
1748 //------------------------------------------------------------------------
1749 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1750 //--------------------------------------------------------------------
1751 // This function calculates the road defined by the cluster density
1752 //--------------------------------------------------------------------
1754 for (Int_t i=0; i<fN; i++) {
1755 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1757 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1759 //------------------------------------------------------------------------
1760 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1761 //--------------------------------------------------------------------
1762 //This function adds a cluster to this layer
1763 //--------------------------------------------------------------------
1764 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1770 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1772 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1773 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1774 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1775 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1776 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1777 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1780 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1781 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1782 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1783 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1787 //------------------------------------------------------------------------
1788 void AliITStrackerMI::AliITSlayer::SortClusters()
1793 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1794 Float_t *z = new Float_t[fN];
1795 Int_t * index = new Int_t[fN];
1797 fMaxSigmaClY=0.; //AD
1798 fMaxSigmaClZ=0.; //AD
1800 for (Int_t i=0;i<fN;i++){
1801 z[i] = fClusters[i]->GetZ();
1802 // save largest errors in y and z for this layer
1803 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1804 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1806 TMath::Sort(fN,z,index,kFALSE);
1807 for (Int_t i=0;i<fN;i++){
1808 clusters[i] = fClusters[index[i]];
1811 for (Int_t i=0;i<fN;i++){
1812 fClusters[i] = clusters[i];
1813 fZ[i] = fClusters[i]->GetZ();
1814 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1815 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1816 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1826 for (Int_t i=0;i<fN;i++){
1827 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1828 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1829 fClusterIndex[i] = i;
1833 fDy5 = (fYB[1]-fYB[0])/5.;
1834 fDy10 = (fYB[1]-fYB[0])/10.;
1835 fDy20 = (fYB[1]-fYB[0])/20.;
1836 for (Int_t i=0;i<6;i++) fN5[i] =0;
1837 for (Int_t i=0;i<11;i++) fN10[i]=0;
1838 for (Int_t i=0;i<21;i++) fN20[i]=0;
1840 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;}
1841 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;}
1842 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;}
1845 for (Int_t i=0;i<fN;i++)
1846 for (Int_t irot=-1;irot<=1;irot++){
1847 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1849 for (Int_t slice=0; slice<6;slice++){
1850 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1851 fClusters5[slice][fN5[slice]] = fClusters[i];
1852 fY5[slice][fN5[slice]] = curY;
1853 fZ5[slice][fN5[slice]] = fZ[i];
1854 fClusterIndex5[slice][fN5[slice]]=i;
1859 for (Int_t slice=0; slice<11;slice++){
1860 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1861 fClusters10[slice][fN10[slice]] = fClusters[i];
1862 fY10[slice][fN10[slice]] = curY;
1863 fZ10[slice][fN10[slice]] = fZ[i];
1864 fClusterIndex10[slice][fN10[slice]]=i;
1869 for (Int_t slice=0; slice<21;slice++){
1870 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1871 fClusters20[slice][fN20[slice]] = fClusters[i];
1872 fY20[slice][fN20[slice]] = curY;
1873 fZ20[slice][fN20[slice]] = fZ[i];
1874 fClusterIndex20[slice][fN20[slice]]=i;
1881 // consistency check
1883 for (Int_t i=0;i<fN-1;i++){
1889 for (Int_t slice=0;slice<21;slice++)
1890 for (Int_t i=0;i<fN20[slice]-1;i++){
1891 if (fZ20[slice][i]>fZ20[slice][i+1]){
1898 //------------------------------------------------------------------------
1899 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1900 //--------------------------------------------------------------------
1901 // This function returns the index of the nearest cluster
1902 //--------------------------------------------------------------------
1905 if (fCurrentSlice<0) {
1914 if (ncl==0) return 0;
1915 Int_t b=0, e=ncl-1, m=(b+e)/2;
1916 for (; b<e; m=(b+e)/2) {
1917 // if (z > fClusters[m]->GetZ()) b=m+1;
1918 if (z > zcl[m]) b=m+1;
1923 //------------------------------------------------------------------------
1924 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 {
1925 //--------------------------------------------------------------------
1926 // This function computes the rectangular road for this track
1927 //--------------------------------------------------------------------
1930 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1931 // take into account the misalignment: propagate track to misaligned detector plane
1932 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1934 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1935 TMath::Sqrt(track->GetSigmaZ2() +
1936 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1937 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1938 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1939 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1940 TMath::Sqrt(track->GetSigmaY2() +
1941 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1942 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1943 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1945 // track at boundary between detectors, enlarge road
1946 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1947 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1948 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1949 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1950 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1951 Float_t tgl = TMath::Abs(track->GetTgl());
1952 if (tgl > 1.) tgl=1.;
1953 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1954 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1955 Float_t snp = TMath::Abs(track->GetSnp());
1956 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1957 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1960 // add to the road a term (up to 2-3 mm) to deal with misalignments
1961 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1962 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1964 Double_t r = fgLayers[ilayer].GetR();
1965 zmin = track->GetZ() - dz;
1966 zmax = track->GetZ() + dz;
1967 ymin = track->GetY() + r*det.GetPhi() - dy;
1968 ymax = track->GetY() + r*det.GetPhi() + dy;
1970 // bring track back to idead detector plane
1971 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1975 //------------------------------------------------------------------------
1976 void AliITStrackerMI::AliITSlayer::
1977 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1978 //--------------------------------------------------------------------
1979 // This function sets the "window"
1980 //--------------------------------------------------------------------
1982 Double_t circle=2*TMath::Pi()*fR;
1988 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1989 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1990 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1991 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1992 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1994 Float_t ymiddle = (fYmin+fYmax)*0.5;
1995 if (ymiddle<fYB[0]) {
1996 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1997 } else if (ymiddle>fYB[1]) {
1998 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
2004 fClustersCs = fClusters;
2005 fClusterIndexCs = fClusterIndex;
2011 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2012 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2013 if (slice<0) slice=0;
2014 if (slice>20) slice=20;
2015 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2017 fCurrentSlice=slice;
2018 fClustersCs = fClusters20[fCurrentSlice];
2019 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2020 fYcs = fY20[fCurrentSlice];
2021 fZcs = fZ20[fCurrentSlice];
2022 fNcs = fN20[fCurrentSlice];
2027 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2028 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2029 if (slice<0) slice=0;
2030 if (slice>10) slice=10;
2031 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2033 fCurrentSlice=slice;
2034 fClustersCs = fClusters10[fCurrentSlice];
2035 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2036 fYcs = fY10[fCurrentSlice];
2037 fZcs = fZ10[fCurrentSlice];
2038 fNcs = fN10[fCurrentSlice];
2043 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2044 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2045 if (slice<0) slice=0;
2046 if (slice>5) slice=5;
2047 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2049 fCurrentSlice=slice;
2050 fClustersCs = fClusters5[fCurrentSlice];
2051 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2052 fYcs = fY5[fCurrentSlice];
2053 fZcs = fZ5[fCurrentSlice];
2054 fNcs = fN5[fCurrentSlice];
2058 fI = FindClusterIndex(fZmin);
2059 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2065 //------------------------------------------------------------------------
2066 Int_t AliITStrackerMI::AliITSlayer::
2067 FindDetectorIndex(Double_t phi, Double_t z) const {
2068 //--------------------------------------------------------------------
2069 //This function finds the detector crossed by the track
2070 //--------------------------------------------------------------------
2072 if (fZOffset<0) // old geometry
2073 dphi = -(phi-fPhiOffset);
2074 else // new geometry
2075 dphi = phi-fPhiOffset;
2078 if (dphi < 0) dphi += 2*TMath::Pi();
2079 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2080 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2081 if (np>=fNladders) np-=fNladders;
2082 if (np<0) np+=fNladders;
2085 Double_t dz=fZOffset-z;
2086 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2087 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2088 if (nz>=fNdetectors || nz<0) {
2089 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2093 // ad hoc correction for 3rd ladder of SDD inner layer,
2094 // which is reversed (rotated by pi around local y)
2095 // this correction is OK only from AliITSv11Hybrid onwards
2096 if (GetR()>12. && GetR()<20.) { // SDD inner
2097 if(np==2) { // 3rd ladder
2098 Double_t posMod252[3];
2099 AliITSgeomTGeo::GetTranslation(252,posMod252);
2100 // check the Z coordinate of Mod 252: if negative
2101 // (old SDD geometry in AliITSv11Hybrid)
2102 // the swap of numeration whould be applied
2103 if(posMod252[2]<0.){
2104 nz = (fNdetectors-1) - nz;
2108 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2111 return np*fNdetectors + nz;
2113 //------------------------------------------------------------------------
2114 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2116 //--------------------------------------------------------------------
2117 // This function returns clusters within the "window"
2118 //--------------------------------------------------------------------
2120 if (fCurrentSlice<0) {
2121 Double_t rpi2 = 2.*fR*TMath::Pi();
2122 for (Int_t i=fI; i<fImax; i++) {
2125 if (fYmax<y) y -= rpi2;
2126 if (fYmin>y) y += rpi2;
2127 if (y<fYmin) continue;
2128 if (y>fYmax) continue;
2130 // skip clusters that are in "extended" road but they
2131 // 3sigma error does not touch the original road
2132 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2133 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2135 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2138 return fClusters[i];
2141 for (Int_t i=fI; i<fImax; i++) {
2142 if (fYcs[i]<fYmin) continue;
2143 if (fYcs[i]>fYmax) continue;
2144 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2145 ci=fClusterIndexCs[i];
2147 return fClustersCs[i];
2152 //------------------------------------------------------------------------
2153 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2155 //--------------------------------------------------------------------
2156 // This function returns the layer thickness at this point (units X0)
2157 //--------------------------------------------------------------------
2159 x0=AliITSRecoParam::GetX0Air();
2160 if (43<fR&&fR<45) { //SSD2
2163 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2164 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2165 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2166 for (Int_t i=0; i<12; i++) {
2167 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2168 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2172 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2173 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2177 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2178 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2181 if (37<fR&&fR<41) { //SSD1
2184 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2185 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2186 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2187 for (Int_t i=0; i<11; i++) {
2188 if (TMath::Abs(z-3.9*i)<0.15) {
2189 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2193 if (TMath::Abs(z+3.9*i)<0.15) {
2194 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2198 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2199 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2202 if (13<fR&&fR<26) { //SDD
2205 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2207 if (TMath::Abs(y-1.80)<0.55) {
2209 for (Int_t j=0; j<20; j++) {
2210 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2211 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2214 if (TMath::Abs(y+1.80)<0.55) {
2216 for (Int_t j=0; j<20; j++) {
2217 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2218 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2222 for (Int_t i=0; i<4; i++) {
2223 if (TMath::Abs(z-7.3*i)<0.60) {
2225 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2228 if (TMath::Abs(z+7.3*i)<0.60) {
2230 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2235 if (6<fR&&fR<8) { //SPD2
2236 Double_t dd=0.0063; x0=21.5;
2238 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2239 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2241 if (3<fR&&fR<5) { //SPD1
2242 Double_t dd=0.0063; x0=21.5;
2244 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2245 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2250 //------------------------------------------------------------------------
2251 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2253 fRmisal(det.fRmisal),
2255 fSinPhi(det.fSinPhi),
2256 fCosPhi(det.fCosPhi),
2262 fNChips(det.fNChips),
2263 fChipIsBad(det.fChipIsBad)
2267 //------------------------------------------------------------------------
2268 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2269 const AliITSDetTypeRec *detTypeRec)
2271 //--------------------------------------------------------------------
2272 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2273 //--------------------------------------------------------------------
2275 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2276 // while in the tracker they start from 0 for each layer
2277 for(Int_t il=0; il<ilayer; il++)
2278 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2281 if (ilayer==0 || ilayer==1) { // ---------- SPD
2283 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2285 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2288 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2292 // Get calibration from AliITSDetTypeRec
2293 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2294 calib->SetModuleIndex(idet);
2295 AliITSCalibration *calibSPDdead = 0;
2296 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2297 if (calib->IsBad() ||
2298 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2301 // printf("lay %d bad %d\n",ilayer,idet);
2304 // Get segmentation from AliITSDetTypeRec
2305 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2307 // Read info about bad chips
2308 fNChips = segm->GetMaximumChipIndex()+1;
2309 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2310 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2311 fChipIsBad = new Bool_t[fNChips];
2312 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2313 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2314 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2315 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2320 //------------------------------------------------------------------------
2321 Double_t AliITStrackerMI::GetEffectiveThickness()
2323 //--------------------------------------------------------------------
2324 // Returns the thickness between the current layer and the vertex (units X0)
2325 //--------------------------------------------------------------------
2328 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2329 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2330 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2334 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2335 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2339 Double_t xn=fgLayers[fI].GetR();
2340 for (Int_t i=0; i<fI; i++) {
2341 Double_t xi=fgLayers[i].GetR();
2342 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2348 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2349 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2352 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2353 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2358 //------------------------------------------------------------------------
2359 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2361 //--------------------------------------------------------------------
2362 // Returns the array of layers between the current layer and the vertex
2363 //--------------------------------------------------------------------
2366 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2367 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2368 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2373 for (int il=fI;il--;) {
2376 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2377 xMS[nl++] = AliITSRecoParam::GetrInsideShield(1);
2380 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2381 xMS[nl++] = AliITSRecoParam::GetrInsideShield(0);
2384 x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2385 xMS[nl++] = fgLayers[il].GetR();
2390 x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2391 xMS[nl++] = AliITSRecoParam::GetrPipe();
2397 //------------------------------------------------------------------------
2398 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2399 //-------------------------------------------------------------------
2400 // This function returns number of clusters within the "window"
2401 //--------------------------------------------------------------------
2403 for (Int_t i=fI; i<fN; i++) {
2404 const AliITSRecPoint *c=fClusters[i];
2405 if (c->GetZ() > fZmax) break;
2406 if (c->IsUsed()) continue;
2407 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2408 Double_t y=fR*det.GetPhi() + c->GetY();
2410 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2411 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2413 if (y<fYmin) continue;
2414 if (y>fYmax) continue;
2419 //------------------------------------------------------------------------
2420 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2421 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2423 //--------------------------------------------------------------------
2424 // This function refits the track "track" at the position "x" using
2425 // the clusters from "clusters"
2426 // If "extra"==kTRUE,
2427 // the clusters from overlapped modules get attached to "track"
2428 // If "planeff"==kTRUE,
2429 // special approach for plane efficiency evaluation is applyed
2430 //--------------------------------------------------------------------
2432 Int_t index[AliITSgeomTGeo::kNLayers];
2434 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2435 Int_t nc=clusters->GetNumberOfClusters();
2436 for (k=0; k<nc; k++) {
2437 Int_t idx=clusters->GetClusterIndex(k);
2438 Int_t ilayer=(idx&0xf0000000)>>28;
2442 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2444 //------------------------------------------------------------------------
2445 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2446 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2448 //--------------------------------------------------------------------
2449 // This function refits the track "track" at the position "x" using
2450 // the clusters from array
2451 // If "extra"==kTRUE,
2452 // the clusters from overlapped modules get attached to "track"
2453 // If "planeff"==kTRUE,
2454 // special approach for plane efficiency evaluation is applyed
2455 //--------------------------------------------------------------------
2456 Int_t index[AliITSgeomTGeo::kNLayers];
2458 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2460 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2461 index[k]=clusters[k];
2464 // special for cosmics and TPC prolonged tracks:
2465 // propagate to the innermost of:
2466 // - innermost layer crossed by the track
2467 // - innermost layer where a cluster was associated to the track
2468 static AliITSRecoParam *repa = NULL;
2470 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2472 repa = AliITSRecoParam::GetHighFluxParam();
2473 AliWarning("Using default AliITSRecoParam class");
2476 Int_t evsp=repa->GetEventSpecie();
2478 if(track->GetESDtrack()) trStatus=track->GetStatus();
2479 Int_t innermostlayer=0;
2480 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2482 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2483 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2484 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2485 index[innermostlayer] >= 0 ) break;
2488 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2491 Int_t modstatus=1; // found
2493 Int_t from, to, step;
2494 if (xx > track->GetX()) {
2495 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2498 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2501 TString dir = (step>0 ? "outward" : "inward");
2503 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2504 AliITSlayer &layer=fgLayers[ilayer];
2505 Double_t r=layer.GetR();
2507 if (step<0 && xx>r) break;
2509 // material between SSD and SDD, SDD and SPD
2510 Double_t hI=ilayer-0.5*step;
2511 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2512 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2513 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2514 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2517 Double_t oldGlobXYZ[3];
2518 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2520 // continue if we are already beyond this layer
2521 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2522 if(step>0 && oldGlobR > r) continue; // going outward
2523 if(step<0 && oldGlobR < r) continue; // going inward
2526 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2528 Int_t idet=layer.FindDetectorIndex(phi,z);
2530 // check if we allow a prolongation without point for large-eta tracks
2531 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2533 modstatus = 4; // out in z
2534 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2535 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2538 // apply correction for material of the current layer
2539 // add time if going outward
2540 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2544 if (idet<0) return kFALSE;
2547 const AliITSdetector &det=layer.GetDetector(idet);
2548 // only for ITS-SA tracks refit
2549 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2551 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2553 track->SetDetectorIndex(idet);
2554 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2556 Double_t dz,zmin,zmax,dy,ymin,ymax;
2558 const AliITSRecPoint *clAcc=0;
2559 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2561 Int_t idx=index[ilayer];
2562 if (idx>=0) { // cluster in this layer
2563 modstatus = 6; // no refit
2564 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2566 if (idet != cl->GetDetectorIndex()) {
2567 idet=cl->GetDetectorIndex();
2568 const AliITSdetector &detc=layer.GetDetector(idet);
2569 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2570 track->SetDetectorIndex(idet);
2571 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2573 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2574 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2578 modstatus = 1; // found
2583 } else { // no cluster in this layer
2585 modstatus = 3; // skipped
2586 // Plane Eff determination:
2587 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2588 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2589 UseTrackForPlaneEff(track,ilayer);
2592 modstatus = 5; // no cls in road
2594 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2595 dz = 0.5*(zmax-zmin);
2596 dy = 0.5*(ymax-ymin);
2597 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2598 if (dead==1) modstatus = 7; // holes in z in SPD
2599 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2604 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2605 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2607 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2610 if (extra && clAcc) { // search for extra clusters in overlapped modules
2611 AliITStrackV2 tmp(*track);
2612 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2613 layer.SelectClusters(zmin,zmax,ymin,ymax);
2615 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2617 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2618 Double_t tolerance=0.1;
2619 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2620 // only clusters in another module! (overlaps)
2621 idetExtra = clExtra->GetDetectorIndex();
2622 if (idet == idetExtra) continue;
2624 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2626 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2627 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2628 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2629 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2631 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2632 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2635 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2636 track->SetExtraModule(ilayer,idetExtra);
2638 } // end search for extra clusters in overlapped modules
2640 // Correct for material of the current layer
2642 // add time if going outward
2643 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2644 track->SetCheckInvariant(kTRUE);
2645 } // end loop on layers
2647 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2651 //------------------------------------------------------------------------
2652 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2655 // calculate normalized chi2
2656 // return NormalizedChi2(track,0);
2659 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2660 // track->fdEdxMismatch=0;
2661 Float_t dedxmismatch =0;
2662 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2664 for (Int_t i = 0;i<6;i++){
2665 if (track->GetClIndex(i)>=0){
2666 Float_t cerry, cerrz;
2667 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2669 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2672 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2673 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2674 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2676 cchi2+=(0.5-ratio)*10.;
2677 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2678 dedxmismatch+=(0.5-ratio)*10.;
2682 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2683 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2684 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2685 if (i<2) chi2+=2*cl->GetDeltaProbability();
2691 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2692 track->SetdEdxMismatch(dedxmismatch);
2696 for (Int_t i = 0;i<4;i++){
2697 if (track->GetClIndex(i)>=0){
2698 Float_t cerry, cerrz;
2699 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2700 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2703 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2704 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2708 for (Int_t i = 4;i<6;i++){
2709 if (track->GetClIndex(i)>=0){
2710 Float_t cerry, cerrz;
2711 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2712 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2715 Float_t cerryb, cerrzb;
2716 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2717 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2720 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2721 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2726 if (track->GetESDtrack()->GetTPCsignal()>85){
2727 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2729 chi2+=(0.5-ratio)*5.;
2732 chi2+=(ratio-2.0)*3;
2736 Double_t match = TMath::Sqrt(track->GetChi22());
2737 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2738 if (!track->GetConstrain()) {
2739 if (track->GetNumberOfClusters()>2) {
2740 match/=track->GetNumberOfClusters()-2.;
2745 if (match<0) match=0;
2747 // penalty factor for missing points (NDeadZone>0), but no penalty
2748 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2749 // or there is a dead from OCDB)
2750 Float_t deadzonefactor = 0.;
2751 if (track->GetNDeadZone()>0.) {
2752 Int_t sumDeadZoneProbability=0;
2753 for(Int_t ilay=0;ilay<6;ilay++) {
2754 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2756 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2757 if(nDeadZoneWithProbNot1>0) {
2758 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2759 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2760 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2762 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2763 deadzonefactor = 3.*(1.1-deadZoneProbability);
2767 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2768 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2769 1./(1.+track->GetNSkipped()));
2770 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2771 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2774 //------------------------------------------------------------------------
2775 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2778 // return matching chi2 between two tracks
2779 Double_t largeChi2=1000.;
2781 AliITStrackMI track3(*track2);
2782 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2784 vec(0,0)=track1->GetY() - track3.GetY();
2785 vec(1,0)=track1->GetZ() - track3.GetZ();
2786 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2787 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2788 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2791 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2792 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2793 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2794 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2795 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2797 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2798 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2799 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2800 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2802 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2803 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2804 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2806 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2807 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2809 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2812 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2813 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2816 //------------------------------------------------------------------------
2817 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2820 // return probability that given point (characterized by z position and error)
2821 // is in SPD dead zone
2822 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2824 Double_t probability = 0.;
2825 Double_t nearestz = 0.,distz=0.;
2826 Int_t nearestzone = -1;
2827 Double_t mindistz = 1000.;
2829 // find closest dead zone
2830 for (Int_t i=0; i<3; i++) {
2831 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2832 if (distz<mindistz) {
2834 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2839 // too far from dead zone
2840 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2843 Double_t zmin, zmax;
2844 if (nearestzone==0) { // dead zone at z = -7
2845 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2846 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2847 } else if (nearestzone==1) { // dead zone at z = 0
2848 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2849 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2850 } else if (nearestzone==2) { // dead zone at z = +7
2851 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2852 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2857 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2859 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2860 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2861 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2864 //------------------------------------------------------------------------
2865 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2868 // calculate normalized chi2
2870 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2872 for (Int_t i = 0;i<6;i++){
2873 if (TMath::Abs(track->GetDy(i))>0){
2874 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2875 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2878 else{chi2[i]=10000;}
2881 TMath::Sort(6,chi2,index,kFALSE);
2882 Float_t max = float(ncl)*fac-1.;
2883 Float_t sumchi=0, sumweight=0;
2884 for (Int_t i=0;i<max+1;i++){
2885 Float_t weight = (i<max)?1.:(max+1.-i);
2886 sumchi+=weight*chi2[index[i]];
2889 Double_t normchi2 = sumchi/sumweight;
2892 //------------------------------------------------------------------------
2893 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2896 // calculate normalized chi2
2897 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2900 for (Int_t i=0;i<6;i++){
2901 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2902 Double_t sy1 = forwardtrack->GetSigmaY(i);
2903 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2904 Double_t sy2 = backtrack->GetSigmaY(i);
2905 Double_t sz2 = backtrack->GetSigmaZ(i);
2906 if (i<2){ sy2=1000.;sz2=1000;}
2908 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2909 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2911 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2912 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2914 res+= nz0*nz0+ny0*ny0;
2917 if (npoints>1) return
2918 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2919 //2*forwardtrack->fNUsed+
2920 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2921 1./(1.+forwardtrack->GetNSkipped()));
2924 //------------------------------------------------------------------------
2925 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2926 //--------------------------------------------------------------------
2927 // Return pointer to a given cluster
2928 //--------------------------------------------------------------------
2929 Int_t l=(index & 0xf0000000) >> 28;
2930 Int_t c=(index & 0x0fffffff) >> 00;
2931 return fgLayers[l].GetWeight(c);
2933 //------------------------------------------------------------------------
2934 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2936 //---------------------------------------------
2937 // register track to the list
2939 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2942 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2943 Int_t index = track->GetClusterIndex(icluster);
2944 Int_t l=(index & 0xf0000000) >> 28;
2945 Int_t c=(index & 0x0fffffff) >> 00;
2946 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2947 for (Int_t itrack=0;itrack<4;itrack++){
2948 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2949 fgLayers[l].SetClusterTracks(itrack,c,id);
2955 //------------------------------------------------------------------------
2956 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2958 //---------------------------------------------
2959 // unregister track from the list
2960 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2961 Int_t index = track->GetClusterIndex(icluster);
2962 Int_t l=(index & 0xf0000000) >> 28;
2963 Int_t c=(index & 0x0fffffff) >> 00;
2964 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2965 for (Int_t itrack=0;itrack<4;itrack++){
2966 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2967 fgLayers[l].SetClusterTracks(itrack,c,-1);
2972 //------------------------------------------------------------------------
2973 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2975 //-------------------------------------------------------------
2976 //get number of shared clusters
2977 //-------------------------------------------------------------
2979 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2980 // mean number of clusters
2981 Float_t *ny = GetNy(id), *nz = GetNz(id);
2984 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2985 Int_t index = track->GetClusterIndex(icluster);
2986 Int_t l=(index & 0xf0000000) >> 28;
2987 Int_t c=(index & 0x0fffffff) >> 00;
2988 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2989 // if (ny[l]<1.e-13){
2990 // printf("problem\n");
2992 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2996 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2997 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2998 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3000 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3003 deltan = (cl->GetNz()-nz[l]);
3005 if (deltan>2.0) continue; // extended - highly probable shared cluster
3006 weight = 2./TMath::Max(3.+deltan,2.);
3008 for (Int_t itrack=0;itrack<4;itrack++){
3009 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3011 clist[l] = (AliITSRecPoint*)GetCluster(index);
3012 track->SetSharedWeight(l,weight);
3018 track->SetNUsed(shared);
3021 //------------------------------------------------------------------------
3022 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3025 // find first shared track
3027 // mean number of clusters
3028 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
3030 for (Int_t i=0;i<6;i++) overlist[i]=-1;
3031 Int_t sharedtrack=100000;
3032 Int_t tracks[24],trackindex=0;
3033 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3035 for (Int_t icluster=0;icluster<6;icluster++){
3036 if (clusterlist[icluster]<0) continue;
3037 Int_t index = clusterlist[icluster];
3038 Int_t l=(index & 0xf0000000) >> 28;
3039 Int_t c=(index & 0x0fffffff) >> 00;
3040 // if (ny[l]<1.e-13){
3041 // printf("problem\n");
3043 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3044 //if (l>3) continue;
3045 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3048 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3049 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3050 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3052 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3055 deltan = (cl->GetNz()-nz[l]);
3057 if (deltan>2.0) continue; // extended - highly probable shared cluster
3059 for (Int_t itrack=3;itrack>=0;itrack--){
3060 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3061 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3062 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
3067 if (trackindex==0) return -1;
3069 sharedtrack = tracks[0];
3071 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3074 Int_t tracks2[24], cluster[24];
3075 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
3078 for (Int_t i=0;i<trackindex;i++){
3079 if (tracks[i]<0) continue;
3080 tracks2[index] = tracks[i];
3082 for (Int_t j=i+1;j<trackindex;j++){
3083 if (tracks[j]<0) continue;
3084 if (tracks[j]==tracks[i]){
3092 for (Int_t i=0;i<index;i++){
3093 if (cluster[index]>max) {
3094 sharedtrack=tracks2[index];
3101 if (sharedtrack>=100000) return -1;
3103 // make list of overlaps
3105 for (Int_t icluster=0;icluster<6;icluster++){
3106 if (clusterlist[icluster]<0) continue;
3107 Int_t index = clusterlist[icluster];
3108 Int_t l=(index & 0xf0000000) >> 28;
3109 Int_t c=(index & 0x0fffffff) >> 00;
3110 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3111 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3113 if (cl->GetNy()>2) continue;
3114 if (cl->GetNz()>2) continue;
3117 if (cl->GetNy()>3) continue;
3118 if (cl->GetNz()>3) continue;
3121 for (Int_t itrack=3;itrack>=0;itrack--){
3122 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3123 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3131 //------------------------------------------------------------------------
3132 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3134 // try to find track hypothesys without conflicts
3135 // with minimal chi2;
3136 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3137 Int_t entries1 = arr1->GetEntriesFast();
3138 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3139 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3140 Int_t entries2 = arr2->GetEntriesFast();
3141 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3143 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3144 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3145 if (track10->Pt()>0.5+track20->Pt()) return track10;
3147 for (Int_t itrack=0;itrack<entries1;itrack++){
3148 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3149 UnRegisterClusterTracks(track,trackID1);
3152 for (Int_t itrack=0;itrack<entries2;itrack++){
3153 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3154 UnRegisterClusterTracks(track,trackID2);
3158 Float_t maxconflicts=6;
3159 Double_t maxchi2 =1000.;
3161 // get weight of hypothesys - using best hypothesy
3164 Int_t list1[6],list2[6];
3165 AliITSRecPoint *clist1[6], *clist2[6] ;
3166 RegisterClusterTracks(track10,trackID1);
3167 RegisterClusterTracks(track20,trackID2);
3168 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3169 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3170 UnRegisterClusterTracks(track10,trackID1);
3171 UnRegisterClusterTracks(track20,trackID2);
3174 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3175 Float_t nerry[6],nerrz[6];
3176 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3177 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3178 for (Int_t i=0;i<6;i++){
3179 if ( (erry1[i]>0) && (erry2[i]>0)) {
3180 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3181 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3183 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3184 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3186 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3187 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3188 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3191 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3192 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3193 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3201 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3202 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3203 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3204 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3206 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3207 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3208 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3210 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3211 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3212 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3215 Double_t sumw = w1+w2;
3219 w1 = (d2+0.5)/(d1+d2+1.);
3220 w2 = (d1+0.5)/(d1+d2+1.);
3222 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3223 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3225 // get pair of "best" hypothesys
3227 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3228 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3230 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3231 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3232 //if (track1->fFakeRatio>0) continue;
3233 RegisterClusterTracks(track1,trackID1);
3234 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3235 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3237 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3238 //if (track2->fFakeRatio>0) continue;
3240 RegisterClusterTracks(track2,trackID2);
3241 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3242 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3243 UnRegisterClusterTracks(track2,trackID2);
3245 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3246 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3247 if (nskipped>0.5) continue;
3249 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3250 if (conflict1+1<cconflict1) continue;
3251 if (conflict2+1<cconflict2) continue;
3255 for (Int_t i=0;i<6;i++){
3257 Float_t c1 =0.; // conflict coeficients
3259 if (clist1[i]&&clist2[i]){
3262 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3265 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3267 c1 = 2./TMath::Max(3.+deltan,2.);
3268 c2 = 2./TMath::Max(3.+deltan,2.);
3274 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3277 deltan = (clist1[i]->GetNz()-nz1[i]);
3279 c1 = 2./TMath::Max(3.+deltan,2.);
3286 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3289 deltan = (clist2[i]->GetNz()-nz2[i]);
3291 c2 = 2./TMath::Max(3.+deltan,2.);
3297 if (TMath::Abs(track1->GetDy(i))>0.) {
3298 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3299 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3300 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3301 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3303 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3306 if (TMath::Abs(track2->GetDy(i))>0.) {
3307 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3308 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3309 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3310 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3313 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3315 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3316 if (chi21>0) sum+=w1;
3317 if (chi22>0) sum+=w2;
3320 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3321 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3322 Double_t normchi2 = 2*conflict+sumchi2/sum;
3323 if ( normchi2 <maxchi2 ){
3326 maxconflicts = conflict;
3330 UnRegisterClusterTracks(track1,trackID1);
3333 // if (maxconflicts<4 && maxchi2<th0){
3334 if (maxchi2<th0*2.){
3335 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3336 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3337 track1->SetChi2MIP(5,maxconflicts);
3338 track1->SetChi2MIP(6,maxchi2);
3339 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3340 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3341 track1->SetChi2MIP(8,index1);
3342 fBestTrackIndex[trackID1] =index1;
3343 UpdateESDtrack(track1, AliESDtrack::kITSin);
3344 original->SetWinner(track1);
3346 else if (track10->GetChi2MIP(0)<th1){
3347 track10->SetChi2MIP(5,maxconflicts);
3348 track10->SetChi2MIP(6,maxchi2);
3349 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3350 UpdateESDtrack(track10,AliESDtrack::kITSin);
3351 original->SetWinner(track10);
3354 for (Int_t itrack=0;itrack<entries1;itrack++){
3355 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3356 UnRegisterClusterTracks(track,trackID1);
3359 for (Int_t itrack=0;itrack<entries2;itrack++){
3360 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3361 UnRegisterClusterTracks(track,trackID2);
3364 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3365 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3366 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3367 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3368 RegisterClusterTracks(track10,trackID1);
3370 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3371 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3372 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3373 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3374 RegisterClusterTracks(track20,trackID2);
3379 //------------------------------------------------------------------------
3380 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3381 //--------------------------------------------------------------------
3382 // This function marks clusters assigned to the track
3383 //--------------------------------------------------------------------
3384 AliTracker::UseClusters(t,from);
3386 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3387 //if (c->GetQ()>2) c->Use();
3388 if (c->GetSigmaZ2()>0.1) c->Use();
3389 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3390 //if (c->GetQ()>2) c->Use();
3391 if (c->GetSigmaZ2()>0.1) c->Use();
3394 //------------------------------------------------------------------------
3395 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3397 //------------------------------------------------------------------
3398 // add track to the list of hypothesys
3399 //------------------------------------------------------------------
3401 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3402 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3404 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3406 array = new TObjArray(10);
3407 fTrackHypothesys.AddAt(array,esdindex);
3409 array->AddLast(track);
3411 //------------------------------------------------------------------------
3412 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3414 //-------------------------------------------------------------------
3415 // compress array of track hypothesys
3416 // keep only maxsize best hypothesys
3417 //-------------------------------------------------------------------
3418 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3419 if (! (fTrackHypothesys.At(esdindex)) ) return;
3420 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3421 Int_t entries = array->GetEntriesFast();
3423 //- find preliminary besttrack as a reference
3424 Float_t minchi2=10000;
3426 AliITStrackMI * besttrack=0;
3428 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3429 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3430 if (!track) continue;
3431 Float_t chi2 = NormalizedChi2(track,0);
3433 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3434 track->SetLabel(tpcLabel);
3436 track->SetFakeRatio(1.);
3437 CookLabel(track,0.); //For comparison only
3439 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3440 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3441 if (track->GetNumberOfClusters()<maxn) continue;
3442 maxn = track->GetNumberOfClusters();
3443 // if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3450 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3451 delete array->RemoveAt(itrack);
3455 if (!besttrack) return;
3458 //take errors of best track as a reference
3459 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3460 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3461 for (Int_t j=0;j<6;j++) {
3462 if (besttrack->GetClIndex(j)>=0){
3463 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3464 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3465 ny[j] = besttrack->GetNy(j);
3466 nz[j] = besttrack->GetNz(j);
3470 // calculate normalized chi2
3472 Float_t * chi2 = new Float_t[entries];
3473 Int_t * index = new Int_t[entries];
3474 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3475 for (Int_t itrack=0;itrack<entries;itrack++){
3476 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3478 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3479 double chi2t = GetNormalizedChi2(track, mode);
3480 track->SetChi2MIP(0,chi2t);
3481 if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3482 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3483 chi2[itrack] = chi2t;
3486 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3487 delete array->RemoveAt(itrack);
3493 TMath::Sort(entries,chi2,index,kFALSE);
3494 besttrack = (AliITStrackMI*)array->At(index[0]);
3495 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3496 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3497 for (Int_t j=0;j<6;j++){
3498 if (besttrack->GetClIndex(j)>=0){
3499 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3500 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3501 ny[j] = besttrack->GetNy(j);
3502 nz[j] = besttrack->GetNz(j);
3507 // calculate one more time with updated normalized errors
3508 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3509 for (Int_t itrack=0;itrack<entries;itrack++){
3510 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3512 double chi2t = GetNormalizedChi2(track, mode);
3513 track->SetChi2MIP(0,chi2t);
3514 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3515 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3516 if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3517 chi2[itrack] = chi2t; //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3520 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3521 delete array->RemoveAt(itrack);
3526 entries = array->GetEntriesFast();
3530 TObjArray * newarray = new TObjArray();
3531 TMath::Sort(entries,chi2,index,kFALSE);
3532 besttrack = (AliITStrackMI*)array->At(index[0]);
3534 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3536 for (Int_t j=0;j<6;j++){
3537 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3538 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3539 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3540 ny[j] = besttrack->GetNy(j);
3541 nz[j] = besttrack->GetNz(j);
3544 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3545 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3546 Float_t minn = besttrack->GetNumberOfClusters()-3;
3548 for (Int_t i=0;i<entries;i++){
3549 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3550 if (!track) continue;
3551 if (accepted>maxcut) break;
3552 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3553 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3554 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3555 delete array->RemoveAt(index[i]);
3559 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3560 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3561 if (!shortbest) accepted++;
3563 newarray->AddLast(array->RemoveAt(index[i]));
3564 for (Int_t j=0;j<6;j++){
3566 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3567 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3568 ny[j] = track->GetNy(j);
3569 nz[j] = track->GetNz(j);
3574 delete array->RemoveAt(index[i]);
3578 delete fTrackHypothesys.RemoveAt(esdindex);
3579 fTrackHypothesys.AddAt(newarray,esdindex);
3583 delete fTrackHypothesys.RemoveAt(esdindex);
3589 //------------------------------------------------------------------------
3590 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3592 //-------------------------------------------------------------
3593 // try to find best hypothesy
3594 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3595 // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3596 //-------------------------------------------------------------
3597 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3598 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3599 if (!array) return 0;
3600 Int_t entries = array->GetEntriesFast();
3601 if (!entries) return 0;
3602 Float_t minchi2 = 100000;
3603 AliITStrackMI * besttrack=0;
3605 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3606 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3607 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3608 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3610 for (Int_t i=0;i<entries;i++){
3611 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3612 if (!track) continue;
3613 Float_t sigmarfi,sigmaz;
3614 GetDCASigma(track,sigmarfi,sigmaz);
3615 track->SetDnorm(0,sigmarfi);
3616 track->SetDnorm(1,sigmaz);
3618 track->SetChi2MIP(1,1000000);
3619 track->SetChi2MIP(2,1000000);
3620 track->SetChi2MIP(3,1000000);
3623 backtrack = new(backtrack) AliITStrackMI(*track);
3624 if (track->GetConstrain()) {
3625 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3626 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3627 if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3628 else {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3630 backtrack->ResetCovariance(10.);
3632 backtrack->ResetCovariance(10.);
3634 backtrack->ResetClusters();
3636 Double_t x = original->GetX();
3637 if (!RefitAt(x,backtrack,track)) continue;
3639 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3640 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3641 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3642 track->SetChi22(GetMatchingChi2(backtrack,original));
3644 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3645 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3646 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3649 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3651 //forward track - without constraint
3652 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3653 forwardtrack->ResetClusters();
3655 if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue; // w/o fwd track MIP03 is meaningless
3656 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3657 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3658 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3660 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3661 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3662 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3663 forwardtrack->SetD(0,track->GetD(0));
3664 forwardtrack->SetD(1,track->GetD(1));
3667 AliITSRecPoint* clist[6];
3668 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3669 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3672 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3673 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3674 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3675 track->SetChi2MIP(3,1000);
3678 Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed(); //RS
3679 if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3680 else chi2 += track->GetNUsed();
3682 for (Int_t ichi=0;ichi<5;ichi++){
3683 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3685 if (chi2 < minchi2){
3686 //besttrack = new AliITStrackMI(*forwardtrack);
3688 besttrack->SetLabel(track->GetLabel());
3689 besttrack->SetFakeRatio(track->GetFakeRatio());
3691 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3692 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3693 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3697 delete forwardtrack;
3699 if (!besttrack) return 0;
3702 for (Int_t i=0;i<entries;i++){
3703 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3705 if (!track) continue;
3707 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3708 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3709 // RS: don't apply this cut when fSelectBestMIP03 is on
3710 || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3712 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3713 delete array->RemoveAt(i);
3723 SortTrackHypothesys(esdindex,checkmax,1);
3725 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3726 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3727 besttrack = (AliITStrackMI*)array->At(0);
3728 if (!besttrack) return 0;
3729 besttrack->SetChi2MIP(8,0);
3730 fBestTrackIndex[esdindex]=0;
3731 entries = array->GetEntriesFast();
3732 AliITStrackMI *longtrack =0;
3734 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3735 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3736 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3737 if (!track->GetConstrain()) continue;
3738 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3739 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3740 if (track->GetChi2MIP(0)>4.) continue;
3741 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3744 //if (longtrack) besttrack=longtrack;
3746 // RS do shared cluster analysis here only if the new sharing analysis is not requested
3747 //RRR if (fFlagFakes) return besttrack;
3750 AliITSRecPoint * clist[6];
3751 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3752 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3753 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3754 RegisterClusterTracks(besttrack,esdindex);
3761 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3762 if (sharedtrack>=0){
3764 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3766 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3772 if (shared>2.5) return 0;
3773 if (shared>1.0) return besttrack;
3775 // Don't sign clusters if not gold track
3777 if (!besttrack->IsGoldPrimary()) return besttrack;
3778 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3780 if (fConstraint[fPass]){
3784 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3785 for (Int_t i=0;i<6;i++){
3786 Int_t index = besttrack->GetClIndex(i);
3787 if (index<0) continue;
3788 Int_t ilayer = (index & 0xf0000000) >> 28;
3789 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3790 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3792 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3793 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3794 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3795 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3796 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3797 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3799 Bool_t cansign = kTRUE;
3800 for (Int_t itrack=0;itrack<entries; itrack++){
3801 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3802 if (!track) continue;
3803 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3804 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3810 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3811 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3812 if (!c->IsUsed()) c->Use();
3818 //------------------------------------------------------------------------
3819 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3822 // get "best" hypothesys
3825 Int_t nentries = itsTracks.GetEntriesFast();
3826 for (Int_t i=0;i<nentries;i++){
3827 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3828 if (!track) continue;
3829 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3830 if (!array) continue;
3831 if (array->GetEntriesFast()<=0) continue;
3833 AliITStrackMI* longtrack=0;
3835 Float_t maxchi2=1000;
3836 for (Int_t j=0;j<array->GetEntriesFast();j++){
3837 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3838 if (!trackHyp) continue;
3839 if (trackHyp->GetGoldV0()) {
3840 longtrack = trackHyp; //gold V0 track taken
3843 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3844 Float_t chi2 = trackHyp->GetChi2MIP(0);
3845 if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3846 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3848 if (fAfterV0){ // ??? RS
3849 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3851 if (chi2 > maxchi2) continue;
3852 minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3853 if (fSelectBestMIP03) minn++; // allow next to longest to win
3860 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3861 if (!longtrack) {longtrack = besttrack;}
3862 else besttrack= longtrack;
3866 AliITSRecPoint * clist[6];
3867 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3869 track->SetNUsed(shared);
3870 track->SetNSkipped(besttrack->GetNSkipped());
3871 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3873 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3877 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3878 //if (sharedtrack==-1) sharedtrack=0;
3879 if (sharedtrack>=0) {
3880 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3883 if (besttrack&&fAfterV0) {
3884 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3885 track->SetWinner(besttrack);
3888 if (fConstraint[fPass]) {
3889 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3890 track->SetWinner(besttrack);
3892 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3893 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3894 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3901 //------------------------------------------------------------------------
3902 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3905 // RS: flag those tracks which are suxpected to have fake clusters
3907 const double kThreshPt = 0.5;
3908 AliRefArray *refArr[6];
3910 for (int i=0;i<6;i++) {
3911 int ncl = fgLayers[i].GetNumberOfClusters();
3912 refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3914 Int_t nentries = itsTracks.GetEntriesFast();
3916 // fill cluster->track associations
3917 for (Int_t itr=0;itr<nentries;itr++){
3918 AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3919 if (!track) continue;
3920 AliITStrackMI* trackITS = track->GetWinner();
3921 if (!trackITS) continue;
3922 for (int il=trackITS->GetNumberOfClusters();il--;) {
3923 int idx = trackITS->GetClusterIndex(il);
3924 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3925 // if (c>fgLayers[l].GetNumberOfClusters()) continue;
3926 refArr[l]->AddReference(c, itr);
3930 const UInt_t kMaxRef = 100;
3931 UInt_t crefs[kMaxRef];
3933 // process tracks with shared clusters
3934 for (int itr=0;itr<nentries;itr++){
3935 AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3936 AliITStrackMI* trackH0 = track0->GetWinner();
3937 if (!trackH0) continue;
3938 AliESDtrack* esd0 = track0->GetESDtrack();
3940 for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3941 int idx = trackH0->GetClusterIndex(il);
3942 Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3943 ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3944 if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3945 esd0->SetITSSharedFlag(l);
3946 for (int ir=ncrefs;ir--;) {
3947 if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3948 AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3949 AliITStrackMI* trackH1 = track1->GetWinner();
3950 AliESDtrack* esd1 = track1->GetESDtrack();
3951 esd1->SetITSSharedFlag(l);
3953 double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3954 if (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3955 else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3957 // select the one with smallest chi2's product
3958 res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3959 res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3961 if (res<0) esd1->SetITSFakeFlag(); // esd0 is winner
3962 else esd0->SetITSFakeFlag(); // esd1 is winner
3969 for (int i=6;i--;) delete refArr[i];
3972 //------------------------------------------------------------------------
3973 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3974 //--------------------------------------------------------------------
3975 //This function "cooks" a track label. If label<0, this track is fake.
3976 //--------------------------------------------------------------------
3979 if (track->GetESDtrack()){
3980 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3981 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3982 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3984 track->SetChi2MIP(9,0);
3986 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3987 Int_t cindex = track->GetClusterIndex(i);
3988 Int_t l=(cindex & 0xf0000000) >> 28;
3989 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3991 for (Int_t ind=0;ind<3;ind++){
3992 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3993 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3995 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3998 Int_t nclusters = track->GetNumberOfClusters();
3999 if (nclusters > 0) //PH Some tracks don't have any cluster
4000 track->SetFakeRatio(double(nwrong)/double(nclusters));
4001 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4002 track->SetLabel(-tpcLabel);
4004 track->SetLabel(tpcLabel);
4006 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4009 //------------------------------------------------------------------------
4010 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4012 // Fill the dE/dx in this track
4014 track->SetChi2MIP(9,0);
4015 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4016 Int_t cindex = track->GetClusterIndex(i);
4017 Int_t l=(cindex & 0xf0000000) >> 28;
4018 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4019 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4021 for (Int_t ind=0;ind<3;ind++){
4022 if (cl->GetLabel(ind)==lab) isWrong=0;
4024 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4028 track->CookdEdx(low,up);
4030 //------------------------------------------------------------------------
4031 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4033 // Create some arrays
4035 if (fCoefficients) delete []fCoefficients;
4036 fCoefficients = new Float_t[ntracks*48];
4037 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4039 //------------------------------------------------------------------------
4040 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4043 // Compute predicted chi2
4045 // Take into account the mis-alignment (bring track to cluster plane)
4046 Double_t xTrOrig=track->GetX();
4047 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4048 Float_t erry,errz,covyz;
4049 Float_t theta = track->GetTgl();
4050 Float_t phi = track->GetSnp();
4051 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4052 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4053 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4054 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4055 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4056 // Bring the track back to detector plane in ideal geometry
4057 // [mis-alignment will be accounted for in UpdateMI()]
4058 if (!track->Propagate(xTrOrig)) return 1000.;
4060 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4061 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4063 chi2+=0.5*TMath::Min(delta/2,2.);
4064 chi2+=2.*cluster->GetDeltaProbability();
4067 track->SetNy(layer,ny);
4068 track->SetNz(layer,nz);
4069 track->SetSigmaY(layer,erry);
4070 track->SetSigmaZ(layer, errz);
4071 track->SetSigmaYZ(layer,covyz);
4072 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4073 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4077 //------------------------------------------------------------------------
4078 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4083 Int_t layer = (index & 0xf0000000) >> 28;
4084 track->SetClIndex(layer, index);
4085 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4086 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4087 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4088 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4092 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
4095 // Take into account the mis-alignment (bring track to cluster plane)
4096 Double_t xTrOrig=track->GetX();
4097 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4098 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4099 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4100 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
4102 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4105 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4106 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4107 c.SetSigmaYZ(track->GetSigmaYZ(layer));
4110 Int_t updated = track->UpdateMI(&c,chi2,index);
4111 // Bring the track back to detector plane in ideal geometry
4112 if (!track->Propagate(xTrOrig)) return 0;
4114 if(!updated) AliDebug(2,"update failed");
4118 //------------------------------------------------------------------------
4119 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4122 //DCA sigmas parameterization
4123 //to be paramterized using external parameters in future
4126 Double_t curv=track->GetC();
4127 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4128 sigmaz = 0.0110+4.37*TMath::Abs(curv);
4130 //------------------------------------------------------------------------
4131 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4134 // Clusters from delta electrons?
4136 Int_t entries = clusterArray->GetEntriesFast();
4137 if (entries<4) return;
4138 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4139 Int_t layer = cluster->GetLayer();
4140 if (layer>1) return;
4142 Int_t ncandidates=0;
4143 Float_t r = (layer>0)? 7:4;
4145 for (Int_t i=0;i<entries;i++){
4146 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4147 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4148 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4149 index[ncandidates] = i; //candidate to belong to delta electron track
4151 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4152 cl0->SetDeltaProbability(1);
4158 for (Int_t i=0;i<ncandidates;i++){
4159 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4160 if (cl0->GetDeltaProbability()>0.8) continue;
4163 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4164 sumy=sumz=sumy2=sumyz=sumw=0.0;
4165 for (Int_t j=0;j<ncandidates;j++){
4167 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4169 Float_t dz = cl0->GetZ()-cl1->GetZ();
4170 Float_t dy = cl0->GetY()-cl1->GetY();
4171 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4172 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4173 y[ncl] = cl1->GetY();
4174 z[ncl] = cl1->GetZ();
4175 sumy+= y[ncl]*weight;
4176 sumz+= z[ncl]*weight;
4177 sumy2+=y[ncl]*y[ncl]*weight;
4178 sumyz+=y[ncl]*z[ncl]*weight;
4183 if (ncl<4) continue;
4184 Float_t det = sumw*sumy2 - sumy*sumy;
4186 if (TMath::Abs(det)>0.01){
4187 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
4188 Float_t k = (sumyz*sumw - sumy*sumz)/det;
4189 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4192 Float_t z0 = sumyz/sumy;
4193 delta = TMath::Abs(cl0->GetZ()-z0);
4196 cl0->SetDeltaProbability(1-20.*delta);
4200 //------------------------------------------------------------------------
4201 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4206 track->UpdateESDtrack(flags);
4207 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4208 if (oldtrack) delete oldtrack;
4209 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4210 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4211 // printf("Problem\n");
4214 //------------------------------------------------------------------------
4215 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4217 // Get nearest upper layer close to the point xr.
4218 // rough approximation
4220 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4221 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4223 for (Int_t i=0;i<6;i++){
4224 if (radius<kRadiuses[i]){
4231 //------------------------------------------------------------------------
4232 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4233 //--------------------------------------------------------------------
4234 // Fill a look-up table with mean material
4235 //--------------------------------------------------------------------
4239 Double_t point1[3],point2[3];
4240 Double_t phi,cosphi,sinphi,z;
4241 // 0-5 layers, 6 pipe, 7-8 shields
4242 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4243 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4245 Int_t ifirst=0,ilast=0;
4246 if(material.Contains("Pipe")) {
4248 } else if(material.Contains("Shields")) {
4250 } else if(material.Contains("Layers")) {
4253 Error("BuildMaterialLUT","Wrong layer name\n");
4256 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4257 Double_t param[5]={0.,0.,0.,0.,0.};
4258 for (Int_t i=0; i<n; i++) {
4259 phi = 2.*TMath::Pi()*gRandom->Rndm();
4260 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4261 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4262 point1[0] = rmin[imat]*cosphi;
4263 point1[1] = rmin[imat]*sinphi;
4265 point2[0] = rmax[imat]*cosphi;
4266 point2[1] = rmax[imat]*sinphi;
4268 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4269 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4271 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4273 fxOverX0Layer[imat] = param[1];
4274 fxTimesRhoLayer[imat] = param[0]*param[4];
4275 } else if(imat==6) {
4276 fxOverX0Pipe = param[1];
4277 fxTimesRhoPipe = param[0]*param[4];
4278 } else if(imat==7) {
4279 fxOverX0Shield[0] = param[1];
4280 fxTimesRhoShield[0] = param[0]*param[4];
4281 } else if(imat==8) {
4282 fxOverX0Shield[1] = param[1];
4283 fxTimesRhoShield[1] = param[0]*param[4];
4287 printf("%s\n",material.Data());
4288 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4289 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4290 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4291 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4292 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4293 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4294 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4295 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4296 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4300 //------------------------------------------------------------------------
4301 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4302 TString direction) {
4303 //-------------------------------------------------------------------
4304 // Propagate beyond beam pipe and correct for material
4305 // (material budget in different ways according to fUseTGeo value)
4306 // Add time if going outward (PropagateTo or PropagateToTGeo)
4307 //-------------------------------------------------------------------
4309 // Define budget mode:
4310 // 0: material from AliITSRecoParam (hard coded)
4311 // 1: material from TGeo in one step (on the fly)
4312 // 2: material from lut
4313 // 3: material from TGeo in one step (same for all hypotheses)
4326 if(fTrackingPhase.Contains("Clusters2Tracks"))
4327 { mode=3; } else { mode=1; }
4330 if(fTrackingPhase.Contains("Clusters2Tracks"))
4331 { mode=3; } else { mode=2; }
4337 if(fTrackingPhase.Contains("Default")) mode=0;
4339 Int_t index=fCurrentEsdTrack;
4341 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4342 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4344 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4346 Double_t xOverX0,x0,lengthTimesMeanDensity;
4350 xOverX0 = AliITSRecoParam::GetdPipe();
4351 x0 = AliITSRecoParam::GetX0Be();
4352 lengthTimesMeanDensity = xOverX0*x0;
4353 lengthTimesMeanDensity *= dir;
4354 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4357 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4360 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4361 xOverX0 = fxOverX0Pipe;
4362 lengthTimesMeanDensity = fxTimesRhoPipe;
4363 lengthTimesMeanDensity *= dir;
4364 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4367 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4368 if(fxOverX0PipeTrks[index]<0) {
4369 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4370 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4371 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4372 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4373 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4376 xOverX0 = fxOverX0PipeTrks[index];
4377 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4378 lengthTimesMeanDensity *= dir;
4379 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4385 //------------------------------------------------------------------------
4386 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4388 TString direction) {
4389 //-------------------------------------------------------------------
4390 // Propagate beyond SPD or SDD shield and correct for material
4391 // (material budget in different ways according to fUseTGeo value)
4392 // Add time if going outward (PropagateTo or PropagateToTGeo)
4393 //-------------------------------------------------------------------
4395 // Define budget mode:
4396 // 0: material from AliITSRecoParam (hard coded)
4397 // 1: material from TGeo in steps of X cm (on the fly)
4398 // X = AliITSRecoParam::GetStepSizeTGeo()
4399 // 2: material from lut
4400 // 3: material from TGeo in one step (same for all hypotheses)
4413 if(fTrackingPhase.Contains("Clusters2Tracks"))
4414 { mode=3; } else { mode=1; }
4417 if(fTrackingPhase.Contains("Clusters2Tracks"))
4418 { mode=3; } else { mode=2; }
4424 if(fTrackingPhase.Contains("Default")) mode=0;
4426 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4428 Int_t shieldindex=0;
4429 if (shield.Contains("SDD")) { // SDDouter
4430 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4432 } else if (shield.Contains("SPD")) { // SPDouter
4433 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4436 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4440 // do nothing if we are already beyond the shield
4441 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4442 if(dir<0 && rTrack > rToGo) return 1; // going outward
4443 if(dir>0 && rTrack < rToGo) return 1; // going inward
4447 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4449 Int_t index=2*fCurrentEsdTrack+shieldindex;
4451 Double_t xOverX0,x0,lengthTimesMeanDensity;
4456 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4457 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4458 lengthTimesMeanDensity = xOverX0*x0;
4459 lengthTimesMeanDensity *= dir;
4460 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4463 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4464 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4467 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4468 xOverX0 = fxOverX0Shield[shieldindex];
4469 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4470 lengthTimesMeanDensity *= dir;
4471 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4474 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4475 if(fxOverX0ShieldTrks[index]<0) {
4476 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4477 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4478 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4479 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4480 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4483 xOverX0 = fxOverX0ShieldTrks[index];
4484 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4485 lengthTimesMeanDensity *= dir;
4486 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4492 //------------------------------------------------------------------------
4493 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4495 Double_t oldGlobXYZ[3],
4496 TString direction) {
4497 //-------------------------------------------------------------------
4498 // Propagate beyond layer and correct for material
4499 // (material budget in different ways according to fUseTGeo value)
4500 // Add time if going outward (PropagateTo or PropagateToTGeo)
4501 //-------------------------------------------------------------------
4503 // Define budget mode:
4504 // 0: material from AliITSRecoParam (hard coded)
4505 // 1: material from TGeo in stepsof X cm (on the fly)
4506 // X = AliITSRecoParam::GetStepSizeTGeo()
4507 // 2: material from lut
4508 // 3: material from TGeo in one step (same for all hypotheses)
4521 if(fTrackingPhase.Contains("Clusters2Tracks"))
4522 { mode=3; } else { mode=1; }
4525 if(fTrackingPhase.Contains("Clusters2Tracks"))
4526 { mode=3; } else { mode=2; }
4532 if(fTrackingPhase.Contains("Default")) mode=0;
4534 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4536 Double_t r=fgLayers[layerindex].GetR();
4537 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4539 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4541 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4543 Int_t index=6*fCurrentEsdTrack+layerindex;
4546 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4549 // back before material (no correction)
4551 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4552 if (!t->GetLocalXat(rOld,xOld)) return 0;
4553 if (!t->Propagate(xOld)) return 0;
4557 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4558 lengthTimesMeanDensity = xOverX0*x0;
4559 lengthTimesMeanDensity *= dir;
4560 // Bring the track beyond the material
4561 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4564 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4565 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4568 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4569 xOverX0 = fxOverX0Layer[layerindex];
4570 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4571 lengthTimesMeanDensity *= dir;
4572 // Bring the track beyond the material
4573 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4576 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4577 if(fxOverX0LayerTrks[index]<0) {
4578 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4579 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4580 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4581 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4582 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4583 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4586 xOverX0 = fxOverX0LayerTrks[index];
4587 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4588 lengthTimesMeanDensity *= dir;
4589 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4596 //------------------------------------------------------------------------
4597 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4598 //-----------------------------------------------------------------
4599 // Initialize LUT for storing material for each prolonged track
4600 //-----------------------------------------------------------------
4601 fxOverX0PipeTrks = new Float_t[ntracks];
4602 fxTimesRhoPipeTrks = new Float_t[ntracks];
4603 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4604 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4605 fxOverX0LayerTrks = new Float_t[ntracks*6];
4606 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4608 for(Int_t i=0; i<ntracks; i++) {
4609 fxOverX0PipeTrks[i] = -1.;
4610 fxTimesRhoPipeTrks[i] = -1.;
4612 for(Int_t j=0; j<ntracks*2; j++) {
4613 fxOverX0ShieldTrks[j] = -1.;
4614 fxTimesRhoShieldTrks[j] = -1.;
4616 for(Int_t k=0; k<ntracks*6; k++) {
4617 fxOverX0LayerTrks[k] = -1.;
4618 fxTimesRhoLayerTrks[k] = -1.;
4625 //------------------------------------------------------------------------
4626 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4627 //-----------------------------------------------------------------
4628 // Delete LUT for storing material for each prolonged track
4629 //-----------------------------------------------------------------
4630 if(fxOverX0PipeTrks) {
4631 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4633 if(fxOverX0ShieldTrks) {
4634 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4637 if(fxOverX0LayerTrks) {
4638 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4640 if(fxTimesRhoPipeTrks) {
4641 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4643 if(fxTimesRhoShieldTrks) {
4644 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4646 if(fxTimesRhoLayerTrks) {
4647 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4651 //------------------------------------------------------------------------
4652 void AliITStrackerMI::SetForceSkippingOfLayer() {
4653 //-----------------------------------------------------------------
4654 // Check if we are forced to skip layers
4655 // either we set to skip them in RecoParam
4656 // or they were off during data-taking
4657 //-----------------------------------------------------------------
4659 const AliEventInfo *eventInfo = GetEventInfo();
4661 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4662 fForceSkippingOfLayer[l] = 0;
4664 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4668 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4669 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4671 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4672 } else if(l==2 || l==3) {
4673 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4675 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4681 //------------------------------------------------------------------------
4682 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4683 Int_t ilayer,Int_t idet) const {
4684 //-----------------------------------------------------------------
4685 // This method is used to decide whether to allow a prolongation
4686 // without clusters, because we want to skip the layer.
4687 // In this case the return value is > 0:
4688 // return 1: the user requested to skip a layer
4689 // return 2: track outside z acceptance
4690 //-----------------------------------------------------------------
4692 if (ForceSkippingOfLayer(ilayer)) return 1;
4694 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4696 if (idet<0 && // out in z
4697 ilayer>innerLayCanSkip &&
4698 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4699 // check if track will cross SPD outer layer
4700 Double_t phiAtSPD2,zAtSPD2;
4701 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4702 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4704 return 2; // always allow skipping, changed on 05.11.2009
4709 //------------------------------------------------------------------------
4710 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4711 Int_t ilayer,Int_t idet,
4712 Double_t dz,Double_t dy,
4713 Bool_t noClusters) const {
4714 //-----------------------------------------------------------------
4715 // This method is used to decide whether to allow a prolongation
4716 // without clusters, because there is a dead zone in the road.
4717 // In this case the return value is > 0:
4718 // return 1: dead zone at z=0,+-7cm in SPD
4719 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4720 // return 2: all road is "bad" (dead or noisy) from the OCDB
4721 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4722 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4723 //-----------------------------------------------------------------
4725 // check dead zones at z=0,+-7cm in the SPD
4726 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4727 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4728 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4729 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4730 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4731 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4732 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4733 for (Int_t i=0; i<3; i++)
4734 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4735 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4736 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4740 // check bad zones from OCDB
4741 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4743 if (idet<0) return 0;
4745 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4748 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4749 if (ilayer==0 || ilayer==1) { // ---------- SPD
4751 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4753 detSizeFactorX *= 2.;
4754 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4757 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4758 if (detType==2) segm->SetLayer(ilayer+1);
4759 Float_t detSizeX = detSizeFactorX*segm->Dx();
4760 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4762 // check if the road overlaps with bad chips
4764 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4765 Float_t zlocmin = zloc-dz;
4766 Float_t zlocmax = zloc+dz;
4767 Float_t xlocmin = xloc-dy;
4768 Float_t xlocmax = xloc+dy;
4769 Int_t chipsInRoad[100];
4771 // check if road goes out of detector
4772 Bool_t touchNeighbourDet=kFALSE;
4773 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4774 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4775 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4776 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4777 AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
4779 // check if this detector is bad
4781 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4782 if(!touchNeighbourDet) {
4783 return 2; // all detectors in road are bad
4785 return 3; // at least one is bad
4789 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4790 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4791 if (!nChipsInRoad) return 0;
4793 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4794 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4795 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4796 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4797 if (det.IsChipBad(chipsInRoad[iCh])) {
4805 if(!touchNeighbourDet) {
4806 AliDebug(2,"all bad in road");
4807 return 2; // all chips in road are bad
4809 return 3; // at least a bad chip in road
4814 AliDebug(2,"at least a bad in road");
4815 return 3; // at least a bad chip in road
4819 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4820 || !noClusters) return 0;
4822 // There are no clusters in road: check if there is at least
4823 // a bad SPD pixel or SDD anode or SSD strips on both sides
4825 Int_t idetInITS=idet;
4826 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4828 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4829 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4832 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4836 //------------------------------------------------------------------------
4837 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4838 const AliITStrackMI *track,
4839 Float_t &xloc,Float_t &zloc) const {
4840 //-----------------------------------------------------------------
4841 // Gives position of track in local module ref. frame
4842 //-----------------------------------------------------------------
4847 if(idet<0) return kTRUE; // track out of z acceptance of layer
4849 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4851 Int_t lad = Int_t(idet/ndet) + 1;
4853 Int_t det = idet - (lad-1)*ndet + 1;
4855 Double_t xyzGlob[3],xyzLoc[3];
4857 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4858 // take into account the misalignment: xyz at real detector plane
4859 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4861 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4863 xloc = (Float_t)xyzLoc[0];
4864 zloc = (Float_t)xyzLoc[2];
4868 //------------------------------------------------------------------------
4869 //------------------------------------------------------------------------
4870 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4872 // Method to be optimized further:
4873 // Aim: decide whether a track can be used for PlaneEff evaluation
4874 // the decision is taken based on the track quality at the layer under study
4875 // no information on the clusters on this layer has to be used
4876 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4877 // the cut is done on number of sigmas from the boundaries
4879 // Input: Actual track, layer [0,5] under study
4881 // Return: kTRUE if this is a good track
4883 // it will apply a pre-selection to obtain good quality tracks.
4884 // Here also you will have the possibility to put a control on the
4885 // impact point of the track on the basic block, in order to exclude border regions
4886 // this will be done by calling a proper method of the AliITSPlaneEff class.
4888 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4889 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4891 Int_t index[AliITSgeomTGeo::kNLayers];
4893 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4895 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4896 index[k]=clusters[k];
4900 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4901 AliITSlayer &layer=fgLayers[ilayer];
4902 Double_t r=layer.GetR();
4903 AliITStrackMI tmp(*track);
4905 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4906 Int_t ncl_out=0; Int_t ncl_in=0;
4907 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4908 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4909 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4910 // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4911 if(index[lay]>=0)ncl_out++;
4913 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4914 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4915 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4916 if (index[lay]>=0) ncl_in++;
4918 Int_t ncl=ncl_out+ncl_in;
4919 Bool_t nextout = kFALSE;
4920 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4921 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4922 Bool_t nextin = kFALSE;
4923 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4924 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4925 // maximum number of missing clusters allowed in outermost layers
4926 if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4928 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4929 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4931 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4932 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4933 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4934 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4938 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4939 Int_t idet=layer.FindDetectorIndex(phi,z);
4940 if(idet<0) { AliInfo(Form("cannot find detector"));
4943 // here check if it has good Chi Square.
4945 //propagate to the intersection with the detector plane
4946 const AliITSdetector &det=layer.GetDetector(idet);
4947 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4951 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4952 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4953 if(key>fPlaneEff->Nblock()) return kFALSE;
4954 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4955 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4957 // DEFINITION OF SEARCH ROAD FOR accepting a track
4959 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4960 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4961 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4962 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4963 // done for RecPoints
4965 // exclude tracks at boundary between detectors
4966 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4967 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4968 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4969 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4970 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4971 if ( (locx-dx < blockXmn+boundaryWidth) ||
4972 (locx+dx > blockXmx-boundaryWidth) ||
4973 (locz-dz < blockZmn+boundaryWidth) ||
4974 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4977 //------------------------------------------------------------------------
4978 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4980 // This Method has to be optimized! For the time-being it uses the same criteria
4981 // as those used in the search of extra clusters for overlapping modules.
4983 // Method Purpose: estabilish whether a track has produced a recpoint or not
4984 // in the layer under study (For Plane efficiency)
4986 // inputs: AliITStrackMI* track (pointer to a usable track)
4988 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4989 // traversed by this very track. In details:
4990 // - if a cluster can be associated to the track then call
4991 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4992 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4995 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4996 AliITSlayer &layer=fgLayers[ilayer];
4997 Double_t r=layer.GetR();
4998 AliITStrackMI tmp(*track);
5002 if (!tmp.GetPhiZat(r,phi,z)) return;
5003 Int_t idet=layer.FindDetectorIndex(phi,z);
5005 if(idet<0) { AliInfo(Form("cannot find detector"));
5009 //propagate to the intersection with the detector plane
5010 const AliITSdetector &det=layer.GetDetector(idet);
5011 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5015 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5017 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5018 TMath::Sqrt(tmp.GetSigmaZ2() +
5019 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5020 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5021 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5022 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5023 TMath::Sqrt(tmp.GetSigmaY2() +
5024 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5025 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5026 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5028 // road in global (rphi,z) [i.e. in tracking ref. system]
5029 Double_t zmin = tmp.GetZ() - dz;
5030 Double_t zmax = tmp.GetZ() + dz;
5031 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5032 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5034 // select clusters in road
5035 layer.SelectClusters(zmin,zmax,ymin,ymax);
5037 // Define criteria for track-cluster association
5038 Double_t msz = tmp.GetSigmaZ2() +
5039 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5040 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5041 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5042 Double_t msy = tmp.GetSigmaY2() +
5043 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5044 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5045 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5046 if (tmp.GetConstrain()) {
5047 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5048 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5050 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5051 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5053 msz = 1./msz; // 1/RoadZ^2
5054 msy = 1./msy; // 1/RoadY^2
5057 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5059 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5060 //Double_t tolerance=0.2;
5061 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5062 idetc = cl->GetDetectorIndex();
5063 if(idet!=idetc) continue;
5064 //Int_t ilay = cl->GetLayer();
5066 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5067 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5069 Double_t chi2=tmp.GetPredictedChi2(cl);
5070 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5074 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5076 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5077 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5078 if(key>fPlaneEff->Nblock()) return;
5079 Bool_t found=kFALSE;
5082 while ((cl=layer.GetNextCluster(clidx))!=0) {
5083 idetc = cl->GetDetectorIndex();
5084 if(idet!=idetc) continue;
5085 // here real control to see whether the cluster can be associated to the track.
5086 // cluster not associated to track
5087 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5088 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5089 // calculate track-clusters chi2
5090 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5091 // in particular, the error associated to the cluster
5092 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5094 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5096 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5097 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5098 // track->SetExtraModule(ilayer,idetExtra);
5100 if(!fPlaneEff->UpDatePlaneEff(found,key))
5101 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5102 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5103 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5104 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5105 Int_t cltype[2]={-999,-999};
5108 Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5112 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5113 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5116 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5117 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5118 cltype[0]=layer.GetCluster(ci)->GetNy();
5119 cltype[1]=layer.GetCluster(ci)->GetNz();
5121 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5122 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5123 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5124 // It is computed properly by calling the method
5125 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5127 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5128 //if (tmp.PropagateTo(x,0.,0.)) {
5129 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5130 AliCluster c(*layer.GetCluster(ci));
5131 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5132 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5133 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5134 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5135 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5138 // Compute the angles between the track and the module
5139 // compute the angle "in phi direction", i.e. the angle in the transverse plane
5140 // between the normal to the module and the projection (in the transverse plane) of the
5142 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5143 Float_t tgl = tmp.GetTgl();
5144 Float_t phitr = tmp.GetSnp();
5145 phitr = TMath::ASin(phitr);
5146 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5148 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5149 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5151 alpha = tmp.GetAlpha();
5152 Double_t phiglob = alpha+phitr;
5154 p[0] = TMath::Cos(phiglob);
5155 p[1] = TMath::Sin(phiglob);
5157 TVector3 pvec(p[0],p[1],p[2]);
5158 TVector3 normvec(rot[1],rot[4],rot[7]);
5159 Double_t angle = pvec.Angle(normvec);
5161 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5162 angle *= 180./TMath::Pi();
5165 TVector3 pt(p[0],p[1],0);
5166 TVector3 normt(rot[1],rot[4],0);
5167 Double_t anglet = pt.Angle(normt);
5169 Double_t phiPt = TMath::ATan2(p[1],p[0]);
5170 if(phiPt<0)phiPt+=2.*TMath::Pi();
5171 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5172 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5173 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5174 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
5175 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5176 anglet *= 180./TMath::Pi();
5178 AngleModTrack[2]=(Float_t) angle;
5179 AngleModTrack[0]=(Float_t) anglet;
5180 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5181 AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5182 AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5183 AngleModTrack[1]*=180./TMath::Pi(); // in degree
5185 fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
5190 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5192 // find the MC cluster for the label
5193 return fgLayers[lr].FindClusterForLabel(label,store);
5197 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5199 // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5200 strncpy(patt,"......",6);
5202 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5205 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5206 Int_t cindex = track->GetClusterIndex(i);
5207 Int_t l=(cindex & 0xf0000000) >> 28;
5208 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5210 for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5211 patt[l] = isWrong ? 'f':'c';
5217 //------------------------------------------------------------------------
5218 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5220 //--------------------------------------------------------------------
5223 for (int ic=0;ic<fN;ic++) {
5224 const AliITSRecPoint *cl = GetCluster(ic);
5225 for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5227 if (store) store[nfound] = ic;