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 "AliITSPlaneEff.h"
39 #include "AliITSCalibrationSPD.h"
40 #include "AliITSCalibrationSDD.h"
41 #include "AliITSCalibrationSSD.h"
42 #include "AliCDBEntry.h"
43 #include "AliCDBManager.h"
44 #include "AliAlignObj.h"
45 #include "AliTrackPointArray.h"
46 #include "AliESDVertex.h"
47 #include "AliESDEvent.h"
48 #include "AliESDtrack.h"
51 #include "AliITSChannelStatus.h"
52 #include "AliITSDetTypeRec.h"
53 #include "AliITSRecPoint.h"
54 #include "AliITSgeomTGeo.h"
55 #include "AliITSReconstructor.h"
56 #include "AliITSClusterParam.h"
57 #include "AliITSsegmentation.h"
58 #include "AliITSCalibration.h"
59 #include "AliITSPlaneEffSPD.h"
60 #include "AliITSPlaneEffSDD.h"
61 #include "AliITSPlaneEffSSD.h"
62 #include "AliITSV0Finder.h"
63 #include "AliITStrackerMI.h"
65 ClassImp(AliITStrackerMI)
67 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
69 AliITStrackerMI::AliITStrackerMI():AliTracker(),
79 fLastLayerToTrackTo(0),
82 fTrackingPhase("Default"),
88 fxTimesRhoPipeTrks(0),
89 fxOverX0ShieldTrks(0),
90 fxTimesRhoShieldTrks(0),
92 fxTimesRhoLayerTrks(0),
99 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
100 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
101 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
103 //------------------------------------------------------------------------
104 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
105 fI(AliITSgeomTGeo::GetNLayers()),
114 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
117 fTrackingPhase("Default"),
123 fxTimesRhoPipeTrks(0),
124 fxOverX0ShieldTrks(0),
125 fxTimesRhoShieldTrks(0),
126 fxOverX0LayerTrks(0),
127 fxTimesRhoLayerTrks(0),
129 fITSChannelStatus(0),
132 //--------------------------------------------------------------------
133 //This is the AliITStrackerMI constructor
134 //--------------------------------------------------------------------
136 AliWarning("\"geom\" is actually a dummy argument !");
142 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
143 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
144 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
146 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
147 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
148 Double_t poff=TMath::ATan2(y,x);
150 Double_t r=TMath::Sqrt(x*x + y*y);
152 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
153 r += TMath::Sqrt(x*x + y*y);
154 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
155 r += TMath::Sqrt(x*x + y*y);
156 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
157 r += TMath::Sqrt(x*x + y*y);
160 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
162 for (Int_t j=1; j<nlad+1; j++) {
163 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
164 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
165 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
167 Double_t txyz[3]={0.};
168 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
169 m.LocalToMaster(txyz,xyz);
170 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
171 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
173 if (phi<0) phi+=TMath::TwoPi();
174 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
176 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
177 new(&det) AliITSdetector(r,phi);
178 // compute the real radius (with misalignment)
179 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
181 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
182 mmisal.LocalToMaster(txyz,xyz);
183 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
184 det.SetRmisal(rmisal);
186 } // end loop on detectors
187 } // end loop on ladders
188 } // end loop on layers
191 fI=AliITSgeomTGeo::GetNLayers();
194 fConstraint[0]=1; fConstraint[1]=0;
196 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
197 AliITSReconstructor::GetRecoParam()->GetYVdef(),
198 AliITSReconstructor::GetRecoParam()->GetZVdef()};
199 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
200 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
201 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
202 SetVertex(xyzVtx,ersVtx);
204 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=AliITSRecoParam::GetLayersNotToSkip(i);
205 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
206 for (Int_t i=0;i<100000;i++){
207 fBestTrackIndex[i]=0;
210 // store positions of centre of SPD modules (in z)
212 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
213 fSPDdetzcentre[0] = tr[2];
214 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
215 fSPDdetzcentre[1] = tr[2];
216 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
217 fSPDdetzcentre[2] = tr[2];
218 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
219 fSPDdetzcentre[3] = tr[2];
221 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
222 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
223 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
227 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
228 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
230 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
232 // only for plane efficiency evaluation
233 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
234 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
235 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
236 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane))
237 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
238 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
239 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
240 else fPlaneEff = new AliITSPlaneEffSSD();
241 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
242 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
243 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
246 //------------------------------------------------------------------------
247 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
249 fBestTrack(tracker.fBestTrack),
250 fTrackToFollow(tracker.fTrackToFollow),
251 fTrackHypothesys(tracker.fTrackHypothesys),
252 fBestHypothesys(tracker.fBestHypothesys),
253 fOriginal(tracker.fOriginal),
254 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
255 fPass(tracker.fPass),
256 fAfterV0(tracker.fAfterV0),
257 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
258 fCoefficients(tracker.fCoefficients),
260 fTrackingPhase(tracker.fTrackingPhase),
261 fUseTGeo(tracker.fUseTGeo),
262 fNtracks(tracker.fNtracks),
263 fxOverX0Pipe(tracker.fxOverX0Pipe),
264 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
266 fxTimesRhoPipeTrks(0),
267 fxOverX0ShieldTrks(0),
268 fxTimesRhoShieldTrks(0),
269 fxOverX0LayerTrks(0),
270 fxTimesRhoLayerTrks(0),
271 fDebugStreamer(tracker.fDebugStreamer),
272 fITSChannelStatus(tracker.fITSChannelStatus),
273 fkDetTypeRec(tracker.fkDetTypeRec),
274 fPlaneEff(tracker.fPlaneEff) {
278 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
281 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
282 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
285 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
286 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
289 //------------------------------------------------------------------------
290 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
291 //Assignment operator
292 this->~AliITStrackerMI();
293 new(this) AliITStrackerMI(tracker);
296 //------------------------------------------------------------------------
297 AliITStrackerMI::~AliITStrackerMI()
302 if (fCoefficients) delete [] fCoefficients;
303 DeleteTrksMaterialLUT();
304 if (fDebugStreamer) {
305 //fDebugStreamer->Close();
306 delete fDebugStreamer;
308 if(fITSChannelStatus) delete fITSChannelStatus;
309 if(fPlaneEff) delete fPlaneEff;
311 //------------------------------------------------------------------------
312 void AliITStrackerMI::SetLayersNotToSkip(const Int_t *l) {
313 //--------------------------------------------------------------------
314 //This function set masks of the layers which must be not skipped
315 //--------------------------------------------------------------------
316 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=l[i];
318 //------------------------------------------------------------------------
319 void AliITStrackerMI::ReadBadFromDetTypeRec() {
320 //--------------------------------------------------------------------
321 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
323 //--------------------------------------------------------------------
325 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
327 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
329 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
332 if(fITSChannelStatus) delete fITSChannelStatus;
333 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
335 // ITS detectors and chips
336 Int_t i=0,j=0,k=0,ndet=0;
337 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
338 Int_t nBadDetsPerLayer=0;
339 ndet=AliITSgeomTGeo::GetNDetectors(i);
340 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
341 for (k=1; k<ndet+1; k++) {
342 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
343 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
344 if(det.IsBad()) {nBadDetsPerLayer++;}
345 } // end loop on detectors
346 } // end loop on ladders
347 Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
348 } // end loop on layers
352 //------------------------------------------------------------------------
353 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
354 //--------------------------------------------------------------------
355 //This function loads ITS clusters
356 //--------------------------------------------------------------------
357 TBranch *branch=cTree->GetBranch("ITSRecPoints");
359 Error("LoadClusters"," can't get the branch !\n");
363 static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
364 branch->SetAddress(&clusters);
366 Int_t i=0,j=0,ndet=0;
368 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
369 ndet=fgLayers[i].GetNdetectors();
370 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
371 for (; j<jmax; j++) {
372 if (!cTree->GetEvent(j)) continue;
373 Int_t ncl=clusters->GetEntriesFast();
374 SignDeltas(clusters,GetZ());
377 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
378 detector=c->GetDetectorIndex();
380 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
382 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
385 // add dead zone "virtual" cluster in SPD, if there is a cluster within
386 // zwindow cm from the dead zone
387 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
388 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
389 Int_t lab[4] = {0,0,0,detector};
390 Int_t info[3] = {0,0,i};
391 Float_t q = 0.; // this identifies virtual clusters
392 Float_t hit[5] = {xdead,
394 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
395 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
397 Bool_t local = kTRUE;
398 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
399 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
400 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
401 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
402 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
403 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
404 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
405 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
406 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
407 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
408 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
409 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
410 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
411 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
412 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
413 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
414 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
415 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
416 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
418 } // "virtual" clusters in SPD
422 fgLayers[i].ResetRoad(); //road defined by the cluster density
423 fgLayers[i].SortClusters();
430 //------------------------------------------------------------------------
431 void AliITStrackerMI::UnloadClusters() {
432 //--------------------------------------------------------------------
433 //This function unloads ITS clusters
434 //--------------------------------------------------------------------
435 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
437 //------------------------------------------------------------------------
438 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
439 //--------------------------------------------------------------------
440 // Publishes all pointers to clusters known to the tracker into the
441 // passed object array.
442 // The ownership is not transfered - the caller is not expected to delete
444 //--------------------------------------------------------------------
446 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
447 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
448 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
455 //------------------------------------------------------------------------
456 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
457 //--------------------------------------------------------------------
458 // Correction for the material between the TPC and the ITS
459 //--------------------------------------------------------------------
460 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
461 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
462 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
463 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
464 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
465 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
466 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
467 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
469 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
475 //------------------------------------------------------------------------
476 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
477 //--------------------------------------------------------------------
478 // This functions reconstructs ITS tracks
479 // The clusters must be already loaded !
480 //--------------------------------------------------------------------
483 fTrackingPhase="Clusters2Tracks";
485 TObjArray itsTracks(15000);
487 fEsd = event; // store pointer to the esd
489 // temporary (for cosmics)
490 if(event->GetVertex()) {
491 TString title = event->GetVertex()->GetTitle();
492 if(title.Contains("cosmics")) {
493 Double_t xyz[3]={GetX(),GetY(),GetZ()};
494 Double_t exyz[3]={0.1,0.1,0.1};
500 {/* Read ESD tracks */
501 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
502 Int_t nentr=event->GetNumberOfTracks();
503 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
505 AliESDtrack *esd=event->GetTrack(nentr);
506 // ---- for debugging:
507 //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); }
509 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
510 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
511 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
512 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
515 t=new AliITStrackMI(*esd);
516 } catch (const Char_t *msg) {
517 //Warning("Clusters2Tracks",msg);
521 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
522 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
525 // look at the ESD mass hypothesys !
526 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
528 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
530 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
531 //track - can be V0 according to TPC
533 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
537 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
541 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
545 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
550 t->SetReconstructed(kFALSE);
551 itsTracks.AddLast(t);
552 fOriginal.AddLast(t);
554 } /* End Read ESD tracks */
558 Int_t nentr=itsTracks.GetEntriesFast();
559 fTrackHypothesys.Expand(nentr);
560 fBestHypothesys.Expand(nentr);
561 MakeCoefficients(nentr);
562 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
564 // THE TWO TRACKING PASSES
565 for (fPass=0; fPass<2; fPass++) {
566 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
567 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
568 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
569 if (t==0) continue; //this track has been already tracked
570 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
571 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
572 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
573 if (fConstraint[fPass]) {
574 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
575 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
578 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
579 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
581 ResetTrackToFollow(*t);
584 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
587 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
589 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
590 if (!besttrack) continue;
591 besttrack->SetLabel(tpcLabel);
592 // besttrack->CookdEdx();
594 besttrack->SetFakeRatio(1.);
595 CookLabel(besttrack,0.); //For comparison only
596 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
598 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
600 t->SetReconstructed(kTRUE);
602 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
604 GetBestHypothesysMIP(itsTracks);
605 } // end loop on the two tracking passes
607 if(event->GetNumberOfV0s()>0) UpdateTPCV0(event);
608 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) FindV02(event);
609 //if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
610 //if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
615 Int_t entries = fTrackHypothesys.GetEntriesFast();
616 for (Int_t ientry=0; ientry<entries; ientry++) {
617 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
618 if (array) array->Delete();
619 delete fTrackHypothesys.RemoveAt(ientry);
622 fTrackHypothesys.Delete();
623 fBestHypothesys.Delete();
625 delete [] fCoefficients;
627 DeleteTrksMaterialLUT();
629 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
631 fTrackingPhase="Default";
635 //------------------------------------------------------------------------
636 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
637 //--------------------------------------------------------------------
638 // This functions propagates reconstructed ITS tracks back
639 // The clusters must be loaded !
640 //--------------------------------------------------------------------
641 fTrackingPhase="PropagateBack";
642 Int_t nentr=event->GetNumberOfTracks();
643 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
646 for (Int_t i=0; i<nentr; i++) {
647 AliESDtrack *esd=event->GetTrack(i);
649 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
650 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
654 t=new AliITStrackMI(*esd);
655 } catch (const Char_t *msg) {
656 //Warning("PropagateBack",msg);
660 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
662 ResetTrackToFollow(*t);
664 // propagate to vertex [SR, GSI 17.02.2003]
665 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
666 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
667 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
668 fTrackToFollow.StartTimeIntegral();
669 // from vertex to outside pipe
670 CorrectForPipeMaterial(&fTrackToFollow,"outward");
673 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
674 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
675 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
679 fTrackToFollow.SetLabel(t->GetLabel());
680 //fTrackToFollow.CookdEdx();
681 CookLabel(&fTrackToFollow,0.); //For comparison only
682 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
683 //UseClusters(&fTrackToFollow);
689 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
691 fTrackingPhase="Default";
695 //------------------------------------------------------------------------
696 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
697 //--------------------------------------------------------------------
698 // This functions refits ITS tracks using the
699 // "inward propagated" TPC tracks
700 // The clusters must be loaded !
701 //--------------------------------------------------------------------
702 fTrackingPhase="RefitInward";
704 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) RefitV02(event);
705 //if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
707 Int_t nentr=event->GetNumberOfTracks();
708 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
711 for (Int_t i=0; i<nentr; i++) {
712 AliESDtrack *esd=event->GetTrack(i);
714 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
715 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
716 if (esd->GetStatus()&AliESDtrack::kTPCout)
717 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
721 t=new AliITStrackMI(*esd);
722 } catch (const Char_t *msg) {
723 //Warning("RefitInward",msg);
727 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
728 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
733 ResetTrackToFollow(*t);
734 fTrackToFollow.ResetClusters();
736 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
737 fTrackToFollow.ResetCovariance(10.);
740 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
741 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
742 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
743 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
744 AliDebug(2," refit OK");
745 fTrackToFollow.SetLabel(t->GetLabel());
746 // fTrackToFollow.CookdEdx();
747 CookdEdx(&fTrackToFollow);
749 CookLabel(&fTrackToFollow,0.0); //For comparison only
752 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
753 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
754 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
755 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
756 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
757 Double_t r[3]={0.,0.,0.};
759 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
766 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
768 fTrackingPhase="Default";
772 //------------------------------------------------------------------------
773 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
774 //--------------------------------------------------------------------
775 // Return pointer to a given cluster
776 //--------------------------------------------------------------------
777 Int_t l=(index & 0xf0000000) >> 28;
778 Int_t c=(index & 0x0fffffff) >> 00;
779 return fgLayers[l].GetCluster(c);
781 //------------------------------------------------------------------------
782 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
783 //--------------------------------------------------------------------
784 // Get track space point with index i
785 //--------------------------------------------------------------------
787 Int_t l=(index & 0xf0000000) >> 28;
788 Int_t c=(index & 0x0fffffff) >> 00;
789 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
790 Int_t idet = cl->GetDetectorIndex();
794 cl->GetGlobalXYZ(xyz);
795 cl->GetGlobalCov(cov);
797 p.SetCharge(cl->GetQ());
798 p.SetDriftTime(cl->GetDriftTime());
799 p.SetChargeRatio(cl->GetChargeRatio());
800 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
803 iLayer = AliGeomManager::kSPD1;
806 iLayer = AliGeomManager::kSPD2;
809 iLayer = AliGeomManager::kSDD1;
812 iLayer = AliGeomManager::kSDD2;
815 iLayer = AliGeomManager::kSSD1;
818 iLayer = AliGeomManager::kSSD2;
821 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
824 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
825 p.SetVolumeID((UShort_t)volid);
828 //------------------------------------------------------------------------
829 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
830 AliTrackPoint& p, const AliESDtrack *t) {
831 //--------------------------------------------------------------------
832 // Get track space point with index i
833 // (assign error estimated during the tracking)
834 //--------------------------------------------------------------------
836 Int_t l=(index & 0xf0000000) >> 28;
837 Int_t c=(index & 0x0fffffff) >> 00;
838 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
839 Int_t idet = cl->GetDetectorIndex();
841 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
843 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
845 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
846 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
847 Double_t alpha = t->GetAlpha();
848 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
849 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
850 phi += alpha-det.GetPhi();
851 Float_t tgphi = TMath::Tan(phi);
853 Float_t tgl = t->GetTgl(); // tgl about const along track
854 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
856 Float_t errlocalx,errlocalz;
857 Bool_t addMisalErr=kFALSE;
858 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
862 cl->GetGlobalXYZ(xyz);
863 // cl->GetGlobalCov(cov);
864 Float_t pos[3] = {0.,0.,0.};
865 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
866 tmpcl.GetGlobalCov(cov);
869 p.SetCharge(cl->GetQ());
870 p.SetDriftTime(cl->GetDriftTime());
871 p.SetChargeRatio(cl->GetChargeRatio());
873 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
876 iLayer = AliGeomManager::kSPD1;
879 iLayer = AliGeomManager::kSPD2;
882 iLayer = AliGeomManager::kSDD1;
885 iLayer = AliGeomManager::kSDD2;
888 iLayer = AliGeomManager::kSSD1;
891 iLayer = AliGeomManager::kSSD2;
894 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
897 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
899 p.SetVolumeID((UShort_t)volid);
902 //------------------------------------------------------------------------
903 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
905 //--------------------------------------------------------------------
906 // Follow prolongation tree
907 //--------------------------------------------------------------------
909 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
910 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
913 AliESDtrack * esd = otrack->GetESDtrack();
914 if (esd->GetV0Index(0)>0) {
915 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
916 // mapping of ESD track is different as ITS track in Containers
917 // Need something more stable
918 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
919 for (Int_t i=0;i<3;i++){
920 Int_t index = esd->GetV0Index(i);
922 AliESDv0 * vertex = fEsd->GetV0(index);
923 if (vertex->GetStatus()<0) continue; // rejected V0
925 if (esd->GetSign()>0) {
926 vertex->SetIndex(0,esdindex);
928 vertex->SetIndex(1,esdindex);
932 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
934 bestarray = new TObjArray(5);
935 fBestHypothesys.AddAt(bestarray,esdindex);
939 //setup tree of the prolongations
941 static AliITStrackMI tracks[7][100];
942 AliITStrackMI *currenttrack;
943 static AliITStrackMI currenttrack1;
944 static AliITStrackMI currenttrack2;
945 static AliITStrackMI backuptrack;
947 Int_t nindexes[7][100];
948 Float_t normalizedchi2[100];
949 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
950 otrack->SetNSkipped(0);
951 new (&(tracks[6][0])) AliITStrackMI(*otrack);
953 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
954 Int_t modstatus = 1; // found
958 // follow prolongations
959 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
960 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
963 AliITSlayer &layer=fgLayers[ilayer];
964 Double_t r = layer.GetR();
970 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
972 if (ntracks[ilayer]>=100) break;
973 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
974 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
975 if (ntracks[ilayer]>15+ilayer){
976 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
977 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
980 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
982 // material between SSD and SDD, SDD and SPD
984 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
986 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
990 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
991 Int_t idet=layer.FindDetectorIndex(phi,z);
993 Double_t trackGlobXYZ1[3];
994 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
996 // Get the budget to the primary vertex for the current track being prolonged
997 Double_t budgetToPrimVertex = GetEffectiveThickness();
999 // check if we allow a prolongation without point
1000 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1002 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1003 // propagate to the layer radius
1004 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1005 if(!vtrack->Propagate(xToGo)) continue;
1006 // apply correction for material of the current layer
1007 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1008 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1009 vtrack->SetClIndex(ilayer,-1);
1010 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1011 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1012 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1014 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1019 // track outside layer acceptance in z
1020 if (idet<0) continue;
1022 //propagate to the intersection with the detector plane
1023 const AliITSdetector &det=layer.GetDetector(idet);
1024 new(¤ttrack2) AliITStrackMI(currenttrack1);
1025 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1026 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1027 currenttrack1.SetDetectorIndex(idet);
1028 currenttrack2.SetDetectorIndex(idet);
1029 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1032 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1034 // road in global (rphi,z) [i.e. in tracking ref. system]
1035 Double_t zmin,zmax,ymin,ymax;
1036 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1038 // select clusters in road
1039 layer.SelectClusters(zmin,zmax,ymin,ymax);
1040 //********************
1042 // Define criteria for track-cluster association
1043 Double_t msz = currenttrack1.GetSigmaZ2() +
1044 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1045 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1046 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1047 Double_t msy = currenttrack1.GetSigmaY2() +
1048 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1049 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1050 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1052 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1053 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1055 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1056 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1058 msz = 1./msz; // 1/RoadZ^2
1059 msy = 1./msy; // 1/RoadY^2
1063 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1065 const AliITSRecPoint *cl=0;
1067 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1068 Bool_t deadzoneSPD=kFALSE;
1069 currenttrack = ¤ttrack1;
1071 // check if the road contains a dead zone
1072 Bool_t noClusters = kFALSE;
1073 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1074 if (noClusters) AliDebug(2,"no clusters in road");
1075 Double_t dz=0.5*(zmax-zmin);
1076 Double_t dy=0.5*(ymax-ymin);
1077 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1078 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1079 // create a prolongation without clusters (check also if there are no clusters in the road)
1082 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1083 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1084 updatetrack->SetClIndex(ilayer,-1);
1086 modstatus = 5; // no cls in road
1087 } else if (dead==1) {
1088 modstatus = 7; // holes in z in SPD
1089 } else if (dead==2 || dead==3) {
1090 modstatus = 2; // dead from OCDB
1092 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1093 // apply correction for material of the current layer
1094 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1095 if (constrain) { // apply vertex constrain
1096 updatetrack->SetConstrain(constrain);
1097 Bool_t isPrim = kTRUE;
1098 if (ilayer<4) { // check that it's close to the vertex
1099 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1100 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1101 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1102 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1103 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1105 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1108 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1109 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1110 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1118 // loop over clusters in the road
1119 while ((cl=layer.GetNextCluster(clidx))!=0) {
1120 if (ntracks[ilayer]>95) break; //space for skipped clusters
1121 Bool_t changedet =kFALSE;
1122 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1123 Int_t idetc=cl->GetDetectorIndex();
1125 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1126 // take into account misalignment (bring track to real detector plane)
1127 Double_t xTrOrig = currenttrack->GetX();
1128 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1129 // a first cut on track-cluster distance
1130 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1131 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1132 { // cluster not associated to track
1133 AliDebug(2,"not associated");
1136 // bring track back to ideal detector plane
1137 if (!currenttrack->Propagate(xTrOrig)) continue;
1138 } else { // have to move track to cluster's detector
1139 const AliITSdetector &detc=layer.GetDetector(idetc);
1140 // a first cut on track-cluster distance
1142 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1143 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1144 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1145 continue; // cluster not associated to track
1147 new (&backuptrack) AliITStrackMI(currenttrack2);
1149 currenttrack =¤ttrack2;
1150 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1151 new (currenttrack) AliITStrackMI(backuptrack);
1155 currenttrack->SetDetectorIndex(idetc);
1156 // Get again the budget to the primary vertex
1157 // for the current track being prolonged, if had to change detector
1158 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1161 // calculate track-clusters chi2
1162 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1164 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1165 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1166 if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1167 if (ntracks[ilayer]>=100) continue;
1168 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1169 updatetrack->SetClIndex(ilayer,-1);
1170 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1172 if (cl->GetQ()!=0) { // real cluster
1173 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1174 AliDebug(2,"update failed");
1177 updatetrack->SetSampledEdx(cl->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
1178 modstatus = 1; // found
1179 } else { // virtual cluster in dead zone
1180 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1181 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1182 modstatus = 7; // holes in z in SPD
1186 Float_t xlocnewdet,zlocnewdet;
1187 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1188 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1191 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1193 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1195 // apply correction for material of the current layer
1196 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1198 if (constrain) { // apply vertex constrain
1199 updatetrack->SetConstrain(constrain);
1200 Bool_t isPrim = kTRUE;
1201 if (ilayer<4) { // check that it's close to the vertex
1202 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1203 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1204 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1205 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1206 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1208 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1209 } //apply vertex constrain
1211 } // create new hypothesis
1213 AliDebug(2,"chi2 too large");
1216 } // loop over possible prolongations
1218 // allow one prolongation without clusters
1219 if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1220 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1221 // apply correction for material of the current layer
1222 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1223 vtrack->SetClIndex(ilayer,-1);
1224 modstatus = 3; // skipped
1225 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1226 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1227 vtrack->IncrementNSkipped();
1231 // allow one prolongation without clusters for tracks with |tgl|>1.1
1232 if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) { //big theta - for low flux
1233 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1234 // apply correction for material of the current layer
1235 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1236 vtrack->SetClIndex(ilayer,-1);
1237 modstatus = 3; // skipped
1238 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1239 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1240 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1245 } // loop over tracks in layer ilayer+1
1247 //loop over track candidates for the current layer
1253 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1254 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1255 if (normalizedchi2[itrack] <
1256 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1260 if (constrain) { // constrain
1261 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1263 } else { // !constrain
1264 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1269 // sort tracks by increasing normalized chi2
1270 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1271 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1272 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1273 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1274 } // end loop over layers
1278 // Now select tracks to be kept
1280 Int_t max = constrain ? 20 : 5;
1282 // tracks that reach layer 0 (SPD inner)
1283 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1284 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1285 if (track.GetNumberOfClusters()<2) continue;
1286 if (!constrain && track.GetNormChi2(0) >
1287 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1290 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1293 // tracks that reach layer 1 (SPD outer)
1294 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1295 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1296 if (track.GetNumberOfClusters()<4) continue;
1297 if (!constrain && track.GetNormChi2(1) >
1298 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1299 if (constrain) track.IncrementNSkipped();
1301 track.SetD(0,track.GetD(GetX(),GetY()));
1302 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1303 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1304 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1307 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1310 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1312 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1313 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1314 if (track.GetNumberOfClusters()<3) continue;
1315 if (!constrain && track.GetNormChi2(2) >
1316 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1317 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1319 track.SetD(0,track.GetD(GetX(),GetY()));
1320 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1321 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1322 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1325 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1331 // register best track of each layer - important for V0 finder
1333 for (Int_t ilayer=0;ilayer<5;ilayer++){
1334 if (ntracks[ilayer]==0) continue;
1335 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1336 if (track.GetNumberOfClusters()<1) continue;
1337 CookLabel(&track,0);
1338 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1342 // update TPC V0 information
1344 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1345 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1346 for (Int_t i=0;i<3;i++){
1347 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1348 if (index==0) break;
1349 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1350 if (vertex->GetStatus()<0) continue; // rejected V0
1352 if (otrack->GetSign()>0) {
1353 vertex->SetIndex(0,esdindex);
1356 vertex->SetIndex(1,esdindex);
1358 //find nearest layer with track info
1359 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1360 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1361 Int_t nearest = nearestold;
1362 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1363 if (ntracks[nearest]==0){
1368 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1369 if (nearestold<5&&nearest<5){
1370 Bool_t accept = track.GetNormChi2(nearest)<10;
1372 if (track.GetSign()>0) {
1373 vertex->SetParamP(track);
1374 vertex->Update(fprimvertex);
1375 //vertex->SetIndex(0,track.fESDtrack->GetID());
1376 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1378 vertex->SetParamN(track);
1379 vertex->Update(fprimvertex);
1380 //vertex->SetIndex(1,track.fESDtrack->GetID());
1381 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1383 vertex->SetStatus(vertex->GetStatus()+1);
1385 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1392 //------------------------------------------------------------------------
1393 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1395 //--------------------------------------------------------------------
1398 return fgLayers[layer];
1400 //------------------------------------------------------------------------
1401 AliITStrackerMI::AliITSlayer::AliITSlayer():
1426 //--------------------------------------------------------------------
1427 //default AliITSlayer constructor
1428 //--------------------------------------------------------------------
1429 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1430 fClusterWeight[i]=0;
1431 fClusterTracks[0][i]=-1;
1432 fClusterTracks[1][i]=-1;
1433 fClusterTracks[2][i]=-1;
1434 fClusterTracks[3][i]=-1;
1437 //------------------------------------------------------------------------
1438 AliITStrackerMI::AliITSlayer::
1439 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1464 //--------------------------------------------------------------------
1465 //main AliITSlayer constructor
1466 //--------------------------------------------------------------------
1467 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1468 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1470 //------------------------------------------------------------------------
1471 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1473 fPhiOffset(layer.fPhiOffset),
1474 fNladders(layer.fNladders),
1475 fZOffset(layer.fZOffset),
1476 fNdetectors(layer.fNdetectors),
1477 fDetectors(layer.fDetectors),
1482 fClustersCs(layer.fClustersCs),
1483 fClusterIndexCs(layer.fClusterIndexCs),
1487 fCurrentSlice(layer.fCurrentSlice),
1494 fAccepted(layer.fAccepted),
1498 //------------------------------------------------------------------------
1499 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1500 //--------------------------------------------------------------------
1501 // AliITSlayer destructor
1502 //--------------------------------------------------------------------
1503 delete [] fDetectors;
1504 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1505 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1506 fClusterWeight[i]=0;
1507 fClusterTracks[0][i]=-1;
1508 fClusterTracks[1][i]=-1;
1509 fClusterTracks[2][i]=-1;
1510 fClusterTracks[3][i]=-1;
1513 //------------------------------------------------------------------------
1514 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1515 //--------------------------------------------------------------------
1516 // This function removes loaded clusters
1517 //--------------------------------------------------------------------
1518 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1519 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1520 fClusterWeight[i]=0;
1521 fClusterTracks[0][i]=-1;
1522 fClusterTracks[1][i]=-1;
1523 fClusterTracks[2][i]=-1;
1524 fClusterTracks[3][i]=-1;
1530 //------------------------------------------------------------------------
1531 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1532 //--------------------------------------------------------------------
1533 // This function reset weights of the clusters
1534 //--------------------------------------------------------------------
1535 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1536 fClusterWeight[i]=0;
1537 fClusterTracks[0][i]=-1;
1538 fClusterTracks[1][i]=-1;
1539 fClusterTracks[2][i]=-1;
1540 fClusterTracks[3][i]=-1;
1542 for (Int_t i=0; i<fN;i++) {
1543 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1544 if (cl&&cl->IsUsed()) cl->Use();
1548 //------------------------------------------------------------------------
1549 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1550 //--------------------------------------------------------------------
1551 // This function calculates the road defined by the cluster density
1552 //--------------------------------------------------------------------
1554 for (Int_t i=0; i<fN; i++) {
1555 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1557 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1559 //------------------------------------------------------------------------
1560 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1561 //--------------------------------------------------------------------
1562 //This function adds a cluster to this layer
1563 //--------------------------------------------------------------------
1564 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1565 ::Error("InsertCluster","Too many clusters !\n");
1571 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1572 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1573 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1574 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1575 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1579 //------------------------------------------------------------------------
1580 void AliITStrackerMI::AliITSlayer::SortClusters()
1585 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1586 Float_t *z = new Float_t[fN];
1587 Int_t * index = new Int_t[fN];
1589 for (Int_t i=0;i<fN;i++){
1590 z[i] = fClusters[i]->GetZ();
1592 TMath::Sort(fN,z,index,kFALSE);
1593 for (Int_t i=0;i<fN;i++){
1594 clusters[i] = fClusters[index[i]];
1597 for (Int_t i=0;i<fN;i++){
1598 fClusters[i] = clusters[i];
1599 fZ[i] = fClusters[i]->GetZ();
1600 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1601 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1602 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1612 for (Int_t i=0;i<fN;i++){
1613 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1614 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1615 fClusterIndex[i] = i;
1619 fDy5 = (fYB[1]-fYB[0])/5.;
1620 fDy10 = (fYB[1]-fYB[0])/10.;
1621 fDy20 = (fYB[1]-fYB[0])/20.;
1622 for (Int_t i=0;i<6;i++) fN5[i] =0;
1623 for (Int_t i=0;i<11;i++) fN10[i]=0;
1624 for (Int_t i=0;i<21;i++) fN20[i]=0;
1626 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;}
1627 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;}
1628 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;}
1631 for (Int_t i=0;i<fN;i++)
1632 for (Int_t irot=-1;irot<=1;irot++){
1633 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1635 for (Int_t slice=0; slice<6;slice++){
1636 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1637 fClusters5[slice][fN5[slice]] = fClusters[i];
1638 fY5[slice][fN5[slice]] = curY;
1639 fZ5[slice][fN5[slice]] = fZ[i];
1640 fClusterIndex5[slice][fN5[slice]]=i;
1645 for (Int_t slice=0; slice<11;slice++){
1646 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1647 fClusters10[slice][fN10[slice]] = fClusters[i];
1648 fY10[slice][fN10[slice]] = curY;
1649 fZ10[slice][fN10[slice]] = fZ[i];
1650 fClusterIndex10[slice][fN10[slice]]=i;
1655 for (Int_t slice=0; slice<21;slice++){
1656 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1657 fClusters20[slice][fN20[slice]] = fClusters[i];
1658 fY20[slice][fN20[slice]] = curY;
1659 fZ20[slice][fN20[slice]] = fZ[i];
1660 fClusterIndex20[slice][fN20[slice]]=i;
1667 // consistency check
1669 for (Int_t i=0;i<fN-1;i++){
1675 for (Int_t slice=0;slice<21;slice++)
1676 for (Int_t i=0;i<fN20[slice]-1;i++){
1677 if (fZ20[slice][i]>fZ20[slice][i+1]){
1684 //------------------------------------------------------------------------
1685 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1686 //--------------------------------------------------------------------
1687 // This function returns the index of the nearest cluster
1688 //--------------------------------------------------------------------
1691 if (fCurrentSlice<0) {
1700 if (ncl==0) return 0;
1701 Int_t b=0, e=ncl-1, m=(b+e)/2;
1702 for (; b<e; m=(b+e)/2) {
1703 // if (z > fClusters[m]->GetZ()) b=m+1;
1704 if (z > zcl[m]) b=m+1;
1709 //------------------------------------------------------------------------
1710 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 {
1711 //--------------------------------------------------------------------
1712 // This function computes the rectangular road for this track
1713 //--------------------------------------------------------------------
1716 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1717 // take into account the misalignment: propagate track to misaligned detector plane
1718 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1720 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1721 TMath::Sqrt(track->GetSigmaZ2() +
1722 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1723 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1724 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1725 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1726 TMath::Sqrt(track->GetSigmaY2() +
1727 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1728 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1729 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1731 // track at boundary between detectors, enlarge road
1732 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1733 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1734 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1735 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1736 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1737 Float_t tgl = TMath::Abs(track->GetTgl());
1738 if (tgl > 1.) tgl=1.;
1739 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1740 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1741 Float_t snp = TMath::Abs(track->GetSnp());
1742 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1743 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1746 // add to the road a term (up to 2-3 mm) to deal with misalignments
1747 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1748 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1750 Double_t r = fgLayers[ilayer].GetR();
1751 zmin = track->GetZ() - dz;
1752 zmax = track->GetZ() + dz;
1753 ymin = track->GetY() + r*det.GetPhi() - dy;
1754 ymax = track->GetY() + r*det.GetPhi() + dy;
1756 // bring track back to idead detector plane
1757 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1761 //------------------------------------------------------------------------
1762 void AliITStrackerMI::AliITSlayer::
1763 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1764 //--------------------------------------------------------------------
1765 // This function sets the "window"
1766 //--------------------------------------------------------------------
1768 Double_t circle=2*TMath::Pi()*fR;
1769 fYmin = ymin; fYmax =ymax;
1770 Float_t ymiddle = (fYmin+fYmax)*0.5;
1771 if (ymiddle<fYB[0]) {
1772 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1773 } else if (ymiddle>fYB[1]) {
1774 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1780 fClustersCs = fClusters;
1781 fClusterIndexCs = fClusterIndex;
1787 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1788 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1789 if (slice<0) slice=0;
1790 if (slice>20) slice=20;
1791 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1793 fCurrentSlice=slice;
1794 fClustersCs = fClusters20[fCurrentSlice];
1795 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1796 fYcs = fY20[fCurrentSlice];
1797 fZcs = fZ20[fCurrentSlice];
1798 fNcs = fN20[fCurrentSlice];
1803 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1804 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1805 if (slice<0) slice=0;
1806 if (slice>10) slice=10;
1807 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1809 fCurrentSlice=slice;
1810 fClustersCs = fClusters10[fCurrentSlice];
1811 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1812 fYcs = fY10[fCurrentSlice];
1813 fZcs = fZ10[fCurrentSlice];
1814 fNcs = fN10[fCurrentSlice];
1819 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1820 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1821 if (slice<0) slice=0;
1822 if (slice>5) slice=5;
1823 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1825 fCurrentSlice=slice;
1826 fClustersCs = fClusters5[fCurrentSlice];
1827 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1828 fYcs = fY5[fCurrentSlice];
1829 fZcs = fZ5[fCurrentSlice];
1830 fNcs = fN5[fCurrentSlice];
1834 fI=FindClusterIndex(zmin); fZmax=zmax;
1835 fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1841 //------------------------------------------------------------------------
1842 Int_t AliITStrackerMI::AliITSlayer::
1843 FindDetectorIndex(Double_t phi, Double_t z) const {
1844 //--------------------------------------------------------------------
1845 //This function finds the detector crossed by the track
1846 //--------------------------------------------------------------------
1848 if (fZOffset<0) // old geometry
1849 dphi = -(phi-fPhiOffset);
1850 else // new geometry
1851 dphi = phi-fPhiOffset;
1854 if (dphi < 0) dphi += 2*TMath::Pi();
1855 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1856 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1857 if (np>=fNladders) np-=fNladders;
1858 if (np<0) np+=fNladders;
1861 Double_t dz=fZOffset-z;
1862 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1863 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1864 if (nz>=fNdetectors) return -1;
1865 if (nz<0) return -1;
1867 // ad hoc correction for 3rd ladder of SDD inner layer,
1868 // which is reversed (rotated by pi around local y)
1869 // this correction is OK only from AliITSv11Hybrid onwards
1870 if (GetR()>12. && GetR()<20.) { // SDD inner
1871 if(np==2) { // 3rd ladder
1872 nz = (fNdetectors-1) - nz;
1875 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1878 return np*fNdetectors + nz;
1880 //------------------------------------------------------------------------
1881 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1883 //--------------------------------------------------------------------
1884 // This function returns clusters within the "window"
1885 //--------------------------------------------------------------------
1887 if (fCurrentSlice<0) {
1888 Double_t rpi2 = 2.*fR*TMath::Pi();
1889 for (Int_t i=fI; i<fImax; i++) {
1891 if (fYmax<y) y -= rpi2;
1892 if (fYmin>y) y += rpi2;
1893 if (y<fYmin) continue;
1894 if (y>fYmax) continue;
1895 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1898 return fClusters[i];
1901 for (Int_t i=fI; i<fImax; i++) {
1902 if (fYcs[i]<fYmin) continue;
1903 if (fYcs[i]>fYmax) continue;
1904 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1905 ci=fClusterIndexCs[i];
1907 return fClustersCs[i];
1912 //------------------------------------------------------------------------
1913 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1915 //--------------------------------------------------------------------
1916 // This function returns the layer thickness at this point (units X0)
1917 //--------------------------------------------------------------------
1919 x0=AliITSRecoParam::GetX0Air();
1920 if (43<fR&&fR<45) { //SSD2
1923 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1924 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1925 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1926 for (Int_t i=0; i<12; i++) {
1927 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1928 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1932 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1933 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1937 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1938 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1941 if (37<fR&&fR<41) { //SSD1
1944 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1945 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1946 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1947 for (Int_t i=0; i<11; i++) {
1948 if (TMath::Abs(z-3.9*i)<0.15) {
1949 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1953 if (TMath::Abs(z+3.9*i)<0.15) {
1954 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1958 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1959 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1962 if (13<fR&&fR<26) { //SDD
1965 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1967 if (TMath::Abs(y-1.80)<0.55) {
1969 for (Int_t j=0; j<20; j++) {
1970 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1971 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1974 if (TMath::Abs(y+1.80)<0.55) {
1976 for (Int_t j=0; j<20; j++) {
1977 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1978 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1982 for (Int_t i=0; i<4; i++) {
1983 if (TMath::Abs(z-7.3*i)<0.60) {
1985 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1988 if (TMath::Abs(z+7.3*i)<0.60) {
1990 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1995 if (6<fR&&fR<8) { //SPD2
1996 Double_t dd=0.0063; x0=21.5;
1998 if (TMath::Abs(y-3.08)>0.5) d+=dd;
1999 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2001 if (3<fR&&fR<5) { //SPD1
2002 Double_t dd=0.0063; x0=21.5;
2004 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2005 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2010 //------------------------------------------------------------------------
2011 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2013 fRmisal(det.fRmisal),
2015 fSinPhi(det.fSinPhi),
2016 fCosPhi(det.fCosPhi),
2022 fNChips(det.fNChips),
2023 fChipIsBad(det.fChipIsBad)
2027 //------------------------------------------------------------------------
2028 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2029 const AliITSDetTypeRec *detTypeRec)
2031 //--------------------------------------------------------------------
2032 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2033 //--------------------------------------------------------------------
2035 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2036 // while in the tracker they start from 0 for each layer
2037 for(Int_t il=0; il<ilayer; il++)
2038 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2041 if (ilayer==0 || ilayer==1) { // ---------- SPD
2043 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2045 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2048 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2052 // Get calibration from AliITSDetTypeRec
2053 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2054 calib->SetModuleIndex(idet);
2055 AliITSCalibration *calibSPDdead = 0;
2056 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2057 if (calib->IsBad() ||
2058 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2061 // printf("lay %d bad %d\n",ilayer,idet);
2064 // Get segmentation from AliITSDetTypeRec
2065 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2067 // Read info about bad chips
2068 fNChips = segm->GetMaximumChipIndex()+1;
2069 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2070 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2071 fChipIsBad = new Bool_t[fNChips];
2072 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2073 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2074 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2075 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2080 //------------------------------------------------------------------------
2081 Double_t AliITStrackerMI::GetEffectiveThickness()
2083 //--------------------------------------------------------------------
2084 // Returns the thickness between the current layer and the vertex (units X0)
2085 //--------------------------------------------------------------------
2088 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2089 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2090 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2094 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2095 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2099 Double_t xn=fgLayers[fI].GetR();
2100 for (Int_t i=0; i<fI; i++) {
2101 Double_t xi=fgLayers[i].GetR();
2102 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2108 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2109 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2112 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2113 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2117 //------------------------------------------------------------------------
2118 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2119 //-------------------------------------------------------------------
2120 // This function returns number of clusters within the "window"
2121 //--------------------------------------------------------------------
2123 for (Int_t i=fI; i<fN; i++) {
2124 const AliITSRecPoint *c=fClusters[i];
2125 if (c->GetZ() > fZmax) break;
2126 if (c->IsUsed()) continue;
2127 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2128 Double_t y=fR*det.GetPhi() + c->GetY();
2130 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2131 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2133 if (y<fYmin) continue;
2134 if (y>fYmax) continue;
2139 //------------------------------------------------------------------------
2140 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2141 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2143 //--------------------------------------------------------------------
2144 // This function refits the track "track" at the position "x" using
2145 // the clusters from "clusters"
2146 // If "extra"==kTRUE,
2147 // the clusters from overlapped modules get attached to "track"
2148 // If "planeff"==kTRUE,
2149 // special approach for plane efficiency evaluation is applyed
2150 //--------------------------------------------------------------------
2152 Int_t index[AliITSgeomTGeo::kNLayers];
2154 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2155 Int_t nc=clusters->GetNumberOfClusters();
2156 for (k=0; k<nc; k++) {
2157 Int_t idx=clusters->GetClusterIndex(k);
2158 Int_t ilayer=(idx&0xf0000000)>>28;
2162 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2164 //------------------------------------------------------------------------
2165 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2166 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2168 //--------------------------------------------------------------------
2169 // This function refits the track "track" at the position "x" using
2170 // the clusters from array
2171 // If "extra"==kTRUE,
2172 // the clusters from overlapped modules get attached to "track"
2173 // If "planeff"==kTRUE,
2174 // special approach for plane efficiency evaluation is applyed
2175 //--------------------------------------------------------------------
2176 Int_t index[AliITSgeomTGeo::kNLayers];
2178 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2180 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2181 index[k]=clusters[k];
2184 // special for cosmics: check which the innermost layer crossed
2186 Int_t innermostlayer=5;
2187 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2188 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2189 if(drphi < fgLayers[innermostlayer].GetR()) break;
2191 //printf(" drphi %f innermost %d\n",drphi,innermostlayer);
2193 Int_t modstatus=1; // found
2195 Int_t from, to, step;
2196 if (xx > track->GetX()) {
2197 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2200 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2203 TString dir = (step>0 ? "outward" : "inward");
2205 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2206 AliITSlayer &layer=fgLayers[ilayer];
2207 Double_t r=layer.GetR();
2208 if (step<0 && xx>r) break;
2210 // material between SSD and SDD, SDD and SPD
2211 Double_t hI=ilayer-0.5*step;
2212 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2213 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2214 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2215 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2218 Double_t oldGlobXYZ[3];
2219 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2222 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2224 Int_t idet=layer.FindDetectorIndex(phi,z);
2226 // check if we allow a prolongation without point for large-eta tracks
2227 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2229 modstatus = 4; // out in z
2230 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2231 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2234 // apply correction for material of the current layer
2235 // add time if going outward
2236 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2240 if (idet<0) return kFALSE;
2242 const AliITSdetector &det=layer.GetDetector(idet);
2243 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2245 track->SetDetectorIndex(idet);
2246 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2248 Double_t dz,zmin,zmax,dy,ymin,ymax;
2250 const AliITSRecPoint *clAcc=0;
2251 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2253 Int_t idx=index[ilayer];
2254 if (idx>=0) { // cluster in this layer
2255 modstatus = 6; // no refit
2256 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2258 if (idet != cl->GetDetectorIndex()) {
2259 idet=cl->GetDetectorIndex();
2260 const AliITSdetector &detc=layer.GetDetector(idet);
2261 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2262 track->SetDetectorIndex(idet);
2263 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2265 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2266 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2270 modstatus = 1; // found
2275 } else { // no cluster in this layer
2277 modstatus = 3; // skipped
2278 // Plane Eff determination:
2279 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2280 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2281 UseTrackForPlaneEff(track,ilayer);
2284 modstatus = 5; // no cls in road
2286 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2287 dz = 0.5*(zmax-zmin);
2288 dy = 0.5*(ymax-ymin);
2289 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2290 if (dead==1) modstatus = 7; // holes in z in SPD
2291 if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2296 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2297 track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2299 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2302 if (extra) { // search for extra clusters in overlapped modules
2303 AliITStrackV2 tmp(*track);
2304 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2305 layer.SelectClusters(zmin,zmax,ymin,ymax);
2307 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2309 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2310 Double_t tolerance=0.1;
2311 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2312 // only clusters in another module! (overlaps)
2313 idetExtra = clExtra->GetDetectorIndex();
2314 if (idet == idetExtra) continue;
2316 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2318 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2319 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2320 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2321 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2323 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2324 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2327 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2328 track->SetExtraModule(ilayer,idetExtra);
2330 } // end search for extra clusters in overlapped modules
2332 // Correct for material of the current layer
2334 // add time if going outward
2335 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2337 } // end loop on layers
2339 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2343 //------------------------------------------------------------------------
2344 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2347 // calculate normalized chi2
2348 // return NormalizedChi2(track,0);
2351 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2352 // track->fdEdxMismatch=0;
2353 Float_t dedxmismatch =0;
2354 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2356 for (Int_t i = 0;i<6;i++){
2357 if (track->GetClIndex(i)>=0){
2358 Float_t cerry, cerrz;
2359 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2361 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2364 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2365 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2366 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2368 cchi2+=(0.5-ratio)*10.;
2369 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2370 dedxmismatch+=(0.5-ratio)*10.;
2374 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2375 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2376 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2377 if (i<2) chi2+=2*cl->GetDeltaProbability();
2383 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2384 track->SetdEdxMismatch(dedxmismatch);
2388 for (Int_t i = 0;i<4;i++){
2389 if (track->GetClIndex(i)>=0){
2390 Float_t cerry, cerrz;
2391 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2392 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2395 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2396 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2400 for (Int_t i = 4;i<6;i++){
2401 if (track->GetClIndex(i)>=0){
2402 Float_t cerry, cerrz;
2403 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2404 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2407 Float_t cerryb, cerrzb;
2408 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2409 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2412 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2413 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2418 if (track->GetESDtrack()->GetTPCsignal()>85){
2419 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2421 chi2+=(0.5-ratio)*5.;
2424 chi2+=(ratio-2.0)*3;
2428 Double_t match = TMath::Sqrt(track->GetChi22());
2429 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2430 if (!track->GetConstrain()) {
2431 if (track->GetNumberOfClusters()>2) {
2432 match/=track->GetNumberOfClusters()-2.;
2437 if (match<0) match=0;
2438 Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2439 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2440 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2441 1./(1.+track->GetNSkipped()));
2445 //------------------------------------------------------------------------
2446 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2449 // return matching chi2 between two tracks
2450 Double_t largeChi2=1000.;
2452 AliITStrackMI track3(*track2);
2453 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2455 vec(0,0)=track1->GetY() - track3.GetY();
2456 vec(1,0)=track1->GetZ() - track3.GetZ();
2457 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2458 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2459 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2462 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2463 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2464 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2465 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2466 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2468 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2469 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2470 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2471 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2473 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2474 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2475 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2477 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2478 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2480 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2483 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2484 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2487 //------------------------------------------------------------------------
2488 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2491 // return probability that given point (characterized by z position and error)
2492 // is in SPD dead zone
2494 Double_t probability = 0.;
2495 Double_t absz = TMath::Abs(zpos);
2496 Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2497 0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2498 if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2499 Double_t zmin, zmax;
2500 if (zpos<-6.) { // dead zone at z = -7
2501 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2502 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2503 } else if (zpos>6.) { // dead zone at z = +7
2504 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2505 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2506 } else if (absz<2.) { // dead zone at z = 0
2507 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2508 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2513 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2515 probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2516 TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2519 //------------------------------------------------------------------------
2520 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2523 // calculate normalized chi2
2525 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2527 for (Int_t i = 0;i<6;i++){
2528 if (TMath::Abs(track->GetDy(i))>0){
2529 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2530 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2533 else{chi2[i]=10000;}
2536 TMath::Sort(6,chi2,index,kFALSE);
2537 Float_t max = float(ncl)*fac-1.;
2538 Float_t sumchi=0, sumweight=0;
2539 for (Int_t i=0;i<max+1;i++){
2540 Float_t weight = (i<max)?1.:(max+1.-i);
2541 sumchi+=weight*chi2[index[i]];
2544 Double_t normchi2 = sumchi/sumweight;
2547 //------------------------------------------------------------------------
2548 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2551 // calculate normalized chi2
2552 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2555 for (Int_t i=0;i<6;i++){
2556 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2557 Double_t sy1 = forwardtrack->GetSigmaY(i);
2558 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2559 Double_t sy2 = backtrack->GetSigmaY(i);
2560 Double_t sz2 = backtrack->GetSigmaZ(i);
2561 if (i<2){ sy2=1000.;sz2=1000;}
2563 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2564 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2566 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2567 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2569 res+= nz0*nz0+ny0*ny0;
2572 if (npoints>1) return
2573 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2574 //2*forwardtrack->fNUsed+
2575 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2576 1./(1.+forwardtrack->GetNSkipped()));
2579 //------------------------------------------------------------------------
2580 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2581 //--------------------------------------------------------------------
2582 // Return pointer to a given cluster
2583 //--------------------------------------------------------------------
2584 Int_t l=(index & 0xf0000000) >> 28;
2585 Int_t c=(index & 0x0fffffff) >> 00;
2586 return fgLayers[l].GetWeight(c);
2588 //------------------------------------------------------------------------
2589 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2591 //---------------------------------------------
2592 // register track to the list
2594 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2597 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2598 Int_t index = track->GetClusterIndex(icluster);
2599 Int_t l=(index & 0xf0000000) >> 28;
2600 Int_t c=(index & 0x0fffffff) >> 00;
2601 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2602 for (Int_t itrack=0;itrack<4;itrack++){
2603 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2604 fgLayers[l].SetClusterTracks(itrack,c,id);
2610 //------------------------------------------------------------------------
2611 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2613 //---------------------------------------------
2614 // unregister track from the list
2615 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2616 Int_t index = track->GetClusterIndex(icluster);
2617 Int_t l=(index & 0xf0000000) >> 28;
2618 Int_t c=(index & 0x0fffffff) >> 00;
2619 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2620 for (Int_t itrack=0;itrack<4;itrack++){
2621 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2622 fgLayers[l].SetClusterTracks(itrack,c,-1);
2627 //------------------------------------------------------------------------
2628 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2630 //-------------------------------------------------------------
2631 //get number of shared clusters
2632 //-------------------------------------------------------------
2634 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2635 // mean number of clusters
2636 Float_t *ny = GetNy(id), *nz = GetNz(id);
2639 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2640 Int_t index = track->GetClusterIndex(icluster);
2641 Int_t l=(index & 0xf0000000) >> 28;
2642 Int_t c=(index & 0x0fffffff) >> 00;
2643 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2645 printf("problem\n");
2647 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2651 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2652 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2653 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2655 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2658 deltan = (cl->GetNz()-nz[l]);
2660 if (deltan>2.0) continue; // extended - highly probable shared cluster
2661 weight = 2./TMath::Max(3.+deltan,2.);
2663 for (Int_t itrack=0;itrack<4;itrack++){
2664 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2666 clist[l] = (AliITSRecPoint*)GetCluster(index);
2672 track->SetNUsed(shared);
2675 //------------------------------------------------------------------------
2676 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2679 // find first shared track
2681 // mean number of clusters
2682 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2684 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2685 Int_t sharedtrack=100000;
2686 Int_t tracks[24],trackindex=0;
2687 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2689 for (Int_t icluster=0;icluster<6;icluster++){
2690 if (clusterlist[icluster]<0) continue;
2691 Int_t index = clusterlist[icluster];
2692 Int_t l=(index & 0xf0000000) >> 28;
2693 Int_t c=(index & 0x0fffffff) >> 00;
2695 printf("problem\n");
2697 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2698 //if (l>3) continue;
2699 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2702 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2703 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2704 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2706 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2709 deltan = (cl->GetNz()-nz[l]);
2711 if (deltan>2.0) continue; // extended - highly probable shared cluster
2713 for (Int_t itrack=3;itrack>=0;itrack--){
2714 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2715 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2716 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2721 if (trackindex==0) return -1;
2723 sharedtrack = tracks[0];
2725 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2728 Int_t tracks2[24], cluster[24];
2729 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2732 for (Int_t i=0;i<trackindex;i++){
2733 if (tracks[i]<0) continue;
2734 tracks2[index] = tracks[i];
2736 for (Int_t j=i+1;j<trackindex;j++){
2737 if (tracks[j]<0) continue;
2738 if (tracks[j]==tracks[i]){
2746 for (Int_t i=0;i<index;i++){
2747 if (cluster[index]>max) {
2748 sharedtrack=tracks2[index];
2755 if (sharedtrack>=100000) return -1;
2757 // make list of overlaps
2759 for (Int_t icluster=0;icluster<6;icluster++){
2760 if (clusterlist[icluster]<0) continue;
2761 Int_t index = clusterlist[icluster];
2762 Int_t l=(index & 0xf0000000) >> 28;
2763 Int_t c=(index & 0x0fffffff) >> 00;
2764 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2765 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2767 if (cl->GetNy()>2) continue;
2768 if (cl->GetNz()>2) continue;
2771 if (cl->GetNy()>3) continue;
2772 if (cl->GetNz()>3) continue;
2775 for (Int_t itrack=3;itrack>=0;itrack--){
2776 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2777 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2785 //------------------------------------------------------------------------
2786 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2788 // try to find track hypothesys without conflicts
2789 // with minimal chi2;
2790 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2791 Int_t entries1 = arr1->GetEntriesFast();
2792 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2793 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2794 Int_t entries2 = arr2->GetEntriesFast();
2795 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2797 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2798 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2799 if (track10->Pt()>0.5+track20->Pt()) return track10;
2801 for (Int_t itrack=0;itrack<entries1;itrack++){
2802 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2803 UnRegisterClusterTracks(track,trackID1);
2806 for (Int_t itrack=0;itrack<entries2;itrack++){
2807 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2808 UnRegisterClusterTracks(track,trackID2);
2812 Float_t maxconflicts=6;
2813 Double_t maxchi2 =1000.;
2815 // get weight of hypothesys - using best hypothesy
2818 Int_t list1[6],list2[6];
2819 AliITSRecPoint *clist1[6], *clist2[6] ;
2820 RegisterClusterTracks(track10,trackID1);
2821 RegisterClusterTracks(track20,trackID2);
2822 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2823 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2824 UnRegisterClusterTracks(track10,trackID1);
2825 UnRegisterClusterTracks(track20,trackID2);
2828 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2829 Float_t nerry[6],nerrz[6];
2830 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2831 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2832 for (Int_t i=0;i<6;i++){
2833 if ( (erry1[i]>0) && (erry2[i]>0)) {
2834 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2835 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2837 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2838 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2840 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2841 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2842 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2845 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2846 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2847 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2855 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2856 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2857 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2858 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2860 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2861 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2862 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2864 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2865 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2866 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2869 Double_t sumw = w1+w2;
2873 w1 = (d2+0.5)/(d1+d2+1.);
2874 w2 = (d1+0.5)/(d1+d2+1.);
2876 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2877 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2879 // get pair of "best" hypothesys
2881 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2882 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2884 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2885 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2886 //if (track1->fFakeRatio>0) continue;
2887 RegisterClusterTracks(track1,trackID1);
2888 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2889 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2891 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2892 //if (track2->fFakeRatio>0) continue;
2894 RegisterClusterTracks(track2,trackID2);
2895 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2896 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2897 UnRegisterClusterTracks(track2,trackID2);
2899 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2900 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2901 if (nskipped>0.5) continue;
2903 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2904 if (conflict1+1<cconflict1) continue;
2905 if (conflict2+1<cconflict2) continue;
2909 for (Int_t i=0;i<6;i++){
2911 Float_t c1 =0.; // conflict coeficients
2913 if (clist1[i]&&clist2[i]){
2916 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2919 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2921 c1 = 2./TMath::Max(3.+deltan,2.);
2922 c2 = 2./TMath::Max(3.+deltan,2.);
2928 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2931 deltan = (clist1[i]->GetNz()-nz1[i]);
2933 c1 = 2./TMath::Max(3.+deltan,2.);
2940 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2943 deltan = (clist2[i]->GetNz()-nz2[i]);
2945 c2 = 2./TMath::Max(3.+deltan,2.);
2951 if (TMath::Abs(track1->GetDy(i))>0.) {
2952 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2953 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2954 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2955 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2957 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2960 if (TMath::Abs(track2->GetDy(i))>0.) {
2961 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2962 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2963 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2964 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2967 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2969 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2970 if (chi21>0) sum+=w1;
2971 if (chi22>0) sum+=w2;
2974 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2975 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2976 Double_t normchi2 = 2*conflict+sumchi2/sum;
2977 if ( normchi2 <maxchi2 ){
2980 maxconflicts = conflict;
2984 UnRegisterClusterTracks(track1,trackID1);
2987 // if (maxconflicts<4 && maxchi2<th0){
2988 if (maxchi2<th0*2.){
2989 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
2990 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2991 track1->SetChi2MIP(5,maxconflicts);
2992 track1->SetChi2MIP(6,maxchi2);
2993 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
2994 // track1->UpdateESDtrack(AliESDtrack::kITSin);
2995 track1->SetChi2MIP(8,index1);
2996 fBestTrackIndex[trackID1] =index1;
2997 UpdateESDtrack(track1, AliESDtrack::kITSin);
2999 else if (track10->GetChi2MIP(0)<th1){
3000 track10->SetChi2MIP(5,maxconflicts);
3001 track10->SetChi2MIP(6,maxchi2);
3002 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3003 UpdateESDtrack(track10,AliESDtrack::kITSin);
3006 for (Int_t itrack=0;itrack<entries1;itrack++){
3007 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3008 UnRegisterClusterTracks(track,trackID1);
3011 for (Int_t itrack=0;itrack<entries2;itrack++){
3012 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3013 UnRegisterClusterTracks(track,trackID2);
3016 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3017 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3018 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3019 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3020 RegisterClusterTracks(track10,trackID1);
3022 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3023 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3024 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3025 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3026 RegisterClusterTracks(track20,trackID2);
3031 //------------------------------------------------------------------------
3032 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3033 //--------------------------------------------------------------------
3034 // This function marks clusters assigned to the track
3035 //--------------------------------------------------------------------
3036 AliTracker::UseClusters(t,from);
3038 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3039 //if (c->GetQ()>2) c->Use();
3040 if (c->GetSigmaZ2()>0.1) c->Use();
3041 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3042 //if (c->GetQ()>2) c->Use();
3043 if (c->GetSigmaZ2()>0.1) c->Use();
3046 //------------------------------------------------------------------------
3047 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3049 //------------------------------------------------------------------
3050 // add track to the list of hypothesys
3051 //------------------------------------------------------------------
3053 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3054 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3056 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3058 array = new TObjArray(10);
3059 fTrackHypothesys.AddAt(array,esdindex);
3061 array->AddLast(track);
3063 //------------------------------------------------------------------------
3064 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3066 //-------------------------------------------------------------------
3067 // compress array of track hypothesys
3068 // keep only maxsize best hypothesys
3069 //-------------------------------------------------------------------
3070 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3071 if (! (fTrackHypothesys.At(esdindex)) ) return;
3072 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3073 Int_t entries = array->GetEntriesFast();
3075 //- find preliminary besttrack as a reference
3076 Float_t minchi2=10000;
3078 AliITStrackMI * besttrack=0;
3079 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3080 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3081 if (!track) continue;
3082 Float_t chi2 = NormalizedChi2(track,0);
3084 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3085 track->SetLabel(tpcLabel);
3087 track->SetFakeRatio(1.);
3088 CookLabel(track,0.); //For comparison only
3090 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3091 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3092 if (track->GetNumberOfClusters()<maxn) continue;
3093 maxn = track->GetNumberOfClusters();
3100 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3101 delete array->RemoveAt(itrack);
3105 if (!besttrack) return;
3108 //take errors of best track as a reference
3109 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3110 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3111 for (Int_t j=0;j<6;j++) {
3112 if (besttrack->GetClIndex(j)>=0){
3113 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3114 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3115 ny[j] = besttrack->GetNy(j);
3116 nz[j] = besttrack->GetNz(j);
3120 // calculate normalized chi2
3122 Float_t * chi2 = new Float_t[entries];
3123 Int_t * index = new Int_t[entries];
3124 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3125 for (Int_t itrack=0;itrack<entries;itrack++){
3126 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3128 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3129 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3130 chi2[itrack] = track->GetChi2MIP(0);
3132 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3133 delete array->RemoveAt(itrack);
3139 TMath::Sort(entries,chi2,index,kFALSE);
3140 besttrack = (AliITStrackMI*)array->At(index[0]);
3141 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3142 for (Int_t j=0;j<6;j++){
3143 if (besttrack->GetClIndex(j)>=0){
3144 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3145 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3146 ny[j] = besttrack->GetNy(j);
3147 nz[j] = besttrack->GetNz(j);
3152 // calculate one more time with updated normalized errors
3153 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3154 for (Int_t itrack=0;itrack<entries;itrack++){
3155 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3157 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3158 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3159 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3162 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3163 delete array->RemoveAt(itrack);
3168 entries = array->GetEntriesFast();
3172 TObjArray * newarray = new TObjArray();
3173 TMath::Sort(entries,chi2,index,kFALSE);
3174 besttrack = (AliITStrackMI*)array->At(index[0]);
3177 for (Int_t j=0;j<6;j++){
3178 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3179 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3180 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3181 ny[j] = besttrack->GetNy(j);
3182 nz[j] = besttrack->GetNz(j);
3185 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3186 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3187 Float_t minn = besttrack->GetNumberOfClusters()-3;
3189 for (Int_t i=0;i<entries;i++){
3190 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3191 if (!track) continue;
3192 if (accepted>maxcut) break;
3193 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3194 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3195 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3196 delete array->RemoveAt(index[i]);
3200 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3201 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3202 if (!shortbest) accepted++;
3204 newarray->AddLast(array->RemoveAt(index[i]));
3205 for (Int_t j=0;j<6;j++){
3207 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3208 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3209 ny[j] = track->GetNy(j);
3210 nz[j] = track->GetNz(j);
3215 delete array->RemoveAt(index[i]);
3219 delete fTrackHypothesys.RemoveAt(esdindex);
3220 fTrackHypothesys.AddAt(newarray,esdindex);
3224 delete fTrackHypothesys.RemoveAt(esdindex);
3230 //------------------------------------------------------------------------
3231 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3233 //-------------------------------------------------------------
3234 // try to find best hypothesy
3235 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3236 //-------------------------------------------------------------
3237 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3238 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3239 if (!array) return 0;
3240 Int_t entries = array->GetEntriesFast();
3241 if (!entries) return 0;
3242 Float_t minchi2 = 100000;
3243 AliITStrackMI * besttrack=0;
3245 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3246 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3247 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3248 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3250 for (Int_t i=0;i<entries;i++){
3251 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3252 if (!track) continue;
3253 Float_t sigmarfi,sigmaz;
3254 GetDCASigma(track,sigmarfi,sigmaz);
3255 track->SetDnorm(0,sigmarfi);
3256 track->SetDnorm(1,sigmaz);
3258 track->SetChi2MIP(1,1000000);
3259 track->SetChi2MIP(2,1000000);
3260 track->SetChi2MIP(3,1000000);
3263 backtrack = new(backtrack) AliITStrackMI(*track);
3264 if (track->GetConstrain()) {
3265 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3266 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3267 backtrack->ResetCovariance(10.);
3269 backtrack->ResetCovariance(10.);
3271 backtrack->ResetClusters();
3273 Double_t x = original->GetX();
3274 if (!RefitAt(x,backtrack,track)) continue;
3276 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3277 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3278 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3279 track->SetChi22(GetMatchingChi2(backtrack,original));
3281 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3282 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3283 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3286 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3288 //forward track - without constraint
3289 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3290 forwardtrack->ResetClusters();
3292 RefitAt(x,forwardtrack,track);
3293 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3294 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3295 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3297 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3298 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3299 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3300 forwardtrack->SetD(0,track->GetD(0));
3301 forwardtrack->SetD(1,track->GetD(1));
3304 AliITSRecPoint* clist[6];
3305 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3306 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3309 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3310 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3311 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3312 track->SetChi2MIP(3,1000);
3315 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3317 for (Int_t ichi=0;ichi<5;ichi++){
3318 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3320 if (chi2 < minchi2){
3321 //besttrack = new AliITStrackMI(*forwardtrack);
3323 besttrack->SetLabel(track->GetLabel());
3324 besttrack->SetFakeRatio(track->GetFakeRatio());
3326 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3327 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3328 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3332 delete forwardtrack;
3334 for (Int_t i=0;i<entries;i++){
3335 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3337 if (!track) continue;
3339 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3340 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3341 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3342 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3343 delete array->RemoveAt(i);
3353 SortTrackHypothesys(esdindex,checkmax,1);
3355 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3356 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3357 besttrack = (AliITStrackMI*)array->At(0);
3358 if (!besttrack) return 0;
3359 besttrack->SetChi2MIP(8,0);
3360 fBestTrackIndex[esdindex]=0;
3361 entries = array->GetEntriesFast();
3362 AliITStrackMI *longtrack =0;
3364 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3365 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3366 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3367 if (!track->GetConstrain()) continue;
3368 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3369 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3370 if (track->GetChi2MIP(0)>4.) continue;
3371 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3374 //if (longtrack) besttrack=longtrack;
3377 AliITSRecPoint * clist[6];
3378 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3379 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3380 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3381 RegisterClusterTracks(besttrack,esdindex);
3388 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3389 if (sharedtrack>=0){
3391 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3393 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3399 if (shared>2.5) return 0;
3400 if (shared>1.0) return besttrack;
3402 // Don't sign clusters if not gold track
3404 if (!besttrack->IsGoldPrimary()) return besttrack;
3405 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3407 if (fConstraint[fPass]){
3411 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3412 for (Int_t i=0;i<6;i++){
3413 Int_t index = besttrack->GetClIndex(i);
3414 if (index<0) continue;
3415 Int_t ilayer = (index & 0xf0000000) >> 28;
3416 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3417 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3419 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3420 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3421 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3422 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3423 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3424 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3426 Bool_t cansign = kTRUE;
3427 for (Int_t itrack=0;itrack<entries; itrack++){
3428 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3429 if (!track) continue;
3430 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3431 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3437 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3438 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3439 if (!c->IsUsed()) c->Use();
3445 //------------------------------------------------------------------------
3446 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3449 // get "best" hypothesys
3452 Int_t nentries = itsTracks.GetEntriesFast();
3453 for (Int_t i=0;i<nentries;i++){
3454 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3455 if (!track) continue;
3456 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3457 if (!array) continue;
3458 if (array->GetEntriesFast()<=0) continue;
3460 AliITStrackMI* longtrack=0;
3462 Float_t maxchi2=1000;
3463 for (Int_t j=0;j<array->GetEntriesFast();j++){
3464 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3465 if (!trackHyp) continue;
3466 if (trackHyp->GetGoldV0()) {
3467 longtrack = trackHyp; //gold V0 track taken
3470 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3471 Float_t chi2 = trackHyp->GetChi2MIP(0);
3473 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3475 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3477 if (chi2 > maxchi2) continue;
3478 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3485 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3486 if (!longtrack) {longtrack = besttrack;}
3487 else besttrack= longtrack;
3491 AliITSRecPoint * clist[6];
3492 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3494 track->SetNUsed(shared);
3495 track->SetNSkipped(besttrack->GetNSkipped());
3496 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3498 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3502 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3503 //if (sharedtrack==-1) sharedtrack=0;
3504 if (sharedtrack>=0) {
3505 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3508 if (besttrack&&fAfterV0) {
3509 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3511 if (besttrack&&fConstraint[fPass])
3512 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3513 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3514 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3515 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3521 //------------------------------------------------------------------------
3522 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3523 //--------------------------------------------------------------------
3524 //This function "cooks" a track label. If label<0, this track is fake.
3525 //--------------------------------------------------------------------
3528 if ( track->GetESDtrack()) tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3530 track->SetChi2MIP(9,0);
3532 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3533 Int_t cindex = track->GetClusterIndex(i);
3534 Int_t l=(cindex & 0xf0000000) >> 28;
3535 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3537 for (Int_t ind=0;ind<3;ind++){
3539 if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3540 AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3542 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3545 Int_t nclusters = track->GetNumberOfClusters();
3546 if (nclusters > 0) //PH Some tracks don't have any cluster
3547 track->SetFakeRatio(double(nwrong)/double(nclusters));
3549 if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3551 track->SetLabel(tpcLabel);
3553 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3556 //------------------------------------------------------------------------
3557 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3562 //AliITSRecPoint * clist[6];
3563 // Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3566 track->SetChi2MIP(9,0);
3567 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3568 Int_t cindex = track->GetClusterIndex(i);
3569 Int_t l=(cindex & 0xf0000000) >> 28;
3570 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3571 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3573 for (Int_t ind=0;ind<3;ind++){
3574 if (cl->GetLabel(ind)==lab) isWrong=0;
3576 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3578 //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue; //shared track
3579 //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3580 //if (l<4&& !(cl->GetType()==1)) continue;
3581 dedx[accepted]= track->GetSampledEdx(i);
3582 //dedx[accepted]= track->fNormQ[l];
3590 TMath::Sort(accepted,dedx,indexes,kFALSE);
3593 Double_t nl=low*accepted, nu =up*accepted;
3595 Float_t sumweight =0;
3596 for (Int_t i=0; i<accepted; i++) {
3598 if (i<nl+0.1) weight = TMath::Max(1.-(nl-i),0.);
3599 if (i>nu-1) weight = TMath::Max(nu-i,0.);
3600 sumamp+= dedx[indexes[i]]*weight;
3603 track->SetdEdx(sumamp/sumweight);
3605 //------------------------------------------------------------------------
3606 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3608 // Create some arrays
3610 if (fCoefficients) delete []fCoefficients;
3611 fCoefficients = new Float_t[ntracks*48];
3612 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3614 //------------------------------------------------------------------------
3615 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3618 // Compute predicted chi2
3621 Float_t theta = track->GetTgl();
3622 Float_t phi = track->GetSnp();
3623 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3624 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3625 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()));
3626 // Take into account the mis-alignment (bring track to cluster plane)
3627 Double_t xTrOrig=track->GetX();
3628 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3629 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()));
3630 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3631 // Bring the track back to detector plane in ideal geometry
3632 // [mis-alignment will be accounted for in UpdateMI()]
3633 if (!track->Propagate(xTrOrig)) return 1000.;
3635 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3636 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3638 chi2+=0.5*TMath::Min(delta/2,2.);
3639 chi2+=2.*cluster->GetDeltaProbability();
3642 track->SetNy(layer,ny);
3643 track->SetNz(layer,nz);
3644 track->SetSigmaY(layer,erry);
3645 track->SetSigmaZ(layer, errz);
3646 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3647 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3651 //------------------------------------------------------------------------
3652 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3657 Int_t layer = (index & 0xf0000000) >> 28;
3658 track->SetClIndex(layer, index);
3659 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3660 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3661 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3662 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3666 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
3669 // Take into account the mis-alignment (bring track to cluster plane)
3670 Double_t xTrOrig=track->GetX();
3671 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3672 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3673 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3674 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3676 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3680 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3681 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3684 Int_t updated = track->UpdateMI(&c,chi2,index);
3686 // Bring the track back to detector plane in ideal geometry
3687 if (!track->Propagate(xTrOrig)) return 0;
3689 if(!updated) AliDebug(2,"update failed");
3693 //------------------------------------------------------------------------
3694 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3697 //DCA sigmas parameterization
3698 //to be paramterized using external parameters in future
3701 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3702 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3704 //------------------------------------------------------------------------
3705 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3708 // Clusters from delta electrons?
3710 Int_t entries = clusterArray->GetEntriesFast();
3711 if (entries<4) return;
3712 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3713 Int_t layer = cluster->GetLayer();
3714 if (layer>1) return;
3716 Int_t ncandidates=0;
3717 Float_t r = (layer>0)? 7:4;
3719 for (Int_t i=0;i<entries;i++){
3720 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3721 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3722 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3723 index[ncandidates] = i; //candidate to belong to delta electron track
3725 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3726 cl0->SetDeltaProbability(1);
3732 for (Int_t i=0;i<ncandidates;i++){
3733 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3734 if (cl0->GetDeltaProbability()>0.8) continue;
3737 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3738 sumy=sumz=sumy2=sumyz=sumw=0.0;
3739 for (Int_t j=0;j<ncandidates;j++){
3741 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3743 Float_t dz = cl0->GetZ()-cl1->GetZ();
3744 Float_t dy = cl0->GetY()-cl1->GetY();
3745 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3746 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3747 y[ncl] = cl1->GetY();
3748 z[ncl] = cl1->GetZ();
3749 sumy+= y[ncl]*weight;
3750 sumz+= z[ncl]*weight;
3751 sumy2+=y[ncl]*y[ncl]*weight;
3752 sumyz+=y[ncl]*z[ncl]*weight;
3757 if (ncl<4) continue;
3758 Float_t det = sumw*sumy2 - sumy*sumy;
3760 if (TMath::Abs(det)>0.01){
3761 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3762 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3763 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3766 Float_t z0 = sumyz/sumy;
3767 delta = TMath::Abs(cl0->GetZ()-z0);
3770 cl0->SetDeltaProbability(1-20.*delta);
3774 //------------------------------------------------------------------------
3775 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3780 track->UpdateESDtrack(flags);
3781 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3782 if (oldtrack) delete oldtrack;
3783 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3784 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3785 printf("Problem\n");
3788 //------------------------------------------------------------------------
3789 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3791 // Get nearest upper layer close to the point xr.
3792 // rough approximation
3794 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3795 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3797 for (Int_t i=0;i<6;i++){
3798 if (radius<kRadiuses[i]){
3806 //------------------------------------------------------------------------
3807 void AliITStrackerMI::UpdateTPCV0(const AliESDEvent *event){
3809 //try to update, or reject TPC V0s
3811 Int_t nv0s = event->GetNumberOfV0s();
3812 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3814 for (Int_t i=0;i<nv0s;i++){
3815 AliESDv0 * vertex = event->GetV0(i);
3816 Int_t ip = vertex->GetIndex(0);
3817 Int_t im = vertex->GetIndex(1);
3819 TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3820 TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3821 AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3822 AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3826 if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3827 if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3828 if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3833 if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3834 if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3835 if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3838 if (vertex->GetStatus()==-100) continue;
3840 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
3841 Int_t clayer = GetNearestLayer(xrp); //I.B.
3842 vertex->SetNBefore(clayer); //
3843 vertex->SetChi2Before(9*clayer); //
3844 vertex->SetNAfter(6-clayer); //
3845 vertex->SetChi2After(0); //
3847 if (clayer >1 ){ // calculate chi2 before vertex
3848 Float_t chi2p = 0, chi2m=0;
3851 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3852 if (trackp->GetClIndex(ilayer)>=0){
3853 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3854 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3865 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3866 if (trackm->GetClIndex(ilayer)>=0){
3867 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3868 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3877 vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3878 if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex
3881 if (clayer < 5 ){ // calculate chi2 after vertex
3882 Float_t chi2p = 0, chi2m=0;
3884 if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3885 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3886 if (trackp->GetClIndex(ilayer)>=0){
3887 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3888 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3898 if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3899 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3900 if (trackm->GetClIndex(ilayer)>=0){
3901 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3902 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3911 vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3912 if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS
3917 //------------------------------------------------------------------------
3918 void AliITStrackerMI::FindV02(AliESDEvent *event)
3923 // Cuts on DCA - R dependent
3924 // max distance DCA between 2 tracks cut
3925 // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3927 const Float_t kMaxDist0 = 0.1;
3928 const Float_t kMaxDist1 = 0.1;
3929 const Float_t kMaxDist = 1;
3930 const Float_t kMinPointAngle = 0.85;
3931 const Float_t kMinPointAngle2 = 0.99;
3932 const Float_t kMinR = 0.5;
3933 const Float_t kMaxR = 220;
3934 //const Float_t kCausality0Cut = 0.19;
3935 //const Float_t kLikelihood01Cut = 0.25;
3936 //const Float_t kPointAngleCut = 0.9996;
3937 const Float_t kCausality0Cut = 0.19;
3938 const Float_t kLikelihood01Cut = 0.45;
3939 const Float_t kLikelihood1Cut = 0.5;
3940 const Float_t kCombinedCut = 0.55;
3944 TTreeSRedirector &cstream = *fDebugStreamer;
3945 Int_t ntracks = event->GetNumberOfTracks();
3946 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3947 fOriginal.Expand(ntracks);
3948 fTrackHypothesys.Expand(ntracks);
3949 fBestHypothesys.Expand(ntracks);
3951 AliHelix * helixes = new AliHelix[ntracks+2];
3952 TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain
3953 TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain
3954 TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain
3955 Bool_t * forbidden = new Bool_t [ntracks+2];
3956 Int_t *itsmap = new Int_t [ntracks+2];
3957 Float_t *dist = new Float_t[ntracks+2];
3958 Float_t *normdist0 = new Float_t[ntracks+2];
3959 Float_t *normdist1 = new Float_t[ntracks+2];
3960 Float_t *normdist = new Float_t[ntracks+2];
3961 Float_t *norm = new Float_t[ntracks+2];
3962 Float_t *maxr = new Float_t[ntracks+2];
3963 Float_t *minr = new Float_t[ntracks+2];
3964 Float_t *minPointAngle= new Float_t[ntracks+2];
3966 AliV0 *pvertex = new AliV0;
3967 AliITStrackMI * dummy= new AliITStrackMI;
3969 AliITStrackMI trackat0; //temporary track for DCA calculation
3971 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3973 // make ITS - ESD map
3975 for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3976 itsmap[itrack] = -1;
3977 forbidden[itrack] = kFALSE;
3978 maxr[itrack] = kMaxR;
3979 minr[itrack] = kMinR;
3980 minPointAngle[itrack] = kMinPointAngle;
3982 for (Int_t itrack=0;itrack<nitstracks;itrack++){
3983 AliITStrackMI * original = (AliITStrackMI*)(fOriginal.At(itrack));
3984 Int_t esdindex = original->GetESDtrack()->GetID();
3985 itsmap[esdindex] = itrack;
3988 // create ITS tracks from ESD tracks if not done before
3990 for (Int_t itrack=0;itrack<ntracks;itrack++){
3991 if (itsmap[itrack]>=0) continue;
3992 AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3993 //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3994 //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3995 tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP()); //I.B.
3996 if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3997 // tracks which can reach inner part of ITS
3998 // propagate track to outer its volume - with correction for material
3999 CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
4001 itsmap[itrack] = nitstracks;
4002 fOriginal.AddAt(tpctrack,nitstracks);
4006 // fill temporary arrays
4008 for (Int_t itrack=0;itrack<ntracks;itrack++){
4009 AliESDtrack * esdtrack = event->GetTrack(itrack);
4010 Int_t itsindex = itsmap[itrack];
4011 AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
4012 if (!original) continue;
4013 AliITStrackMI *bestConst = 0;
4014 AliITStrackMI *bestLong = 0;
4015 AliITStrackMI *best = 0;
4018 TObjArray * array = (TObjArray*) fTrackHypothesys.At(itsindex);
4019 Int_t hentries = (array==0) ? 0 : array->GetEntriesFast();
4020 // Get best track with vertex constrain
4021 for (Int_t ih=0;ih<hentries;ih++){
4022 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4023 if (!trackh->GetConstrain()) continue;
4024 if (!bestConst) bestConst = trackh;
4025 if (trackh->GetNumberOfClusters()>5.0){
4026 bestConst = trackh; // full track - with minimal chi2
4029 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()) continue;
4033 // Get best long track without vertex constrain and best track without vertex constrain
4034 for (Int_t ih=0;ih<hentries;ih++){
4035 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4036 if (trackh->GetConstrain()) continue;
4037 if (!best) best = trackh;
4038 if (!bestLong) bestLong = trackh;
4039 if (trackh->GetNumberOfClusters()>5.0){
4040 bestLong = trackh; // full track - with minimal chi2
4043 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()) continue;
4048 bestLong = original;
4050 //I.B. trackat0 = *bestLong;
4051 new (&trackat0) AliITStrackMI(*bestLong);
4052 Double_t xx,yy,zz,alpha;
4053 if (!bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz)) continue;
4054 alpha = TMath::ATan2(yy,xx);
4055 trackat0.Propagate(alpha,0); //PH The check on the return value is temporarily disabled (bug 45751)
4056 // calculate normalized distances to the vertex
4058 Float_t ptfac = (1.+100.*TMath::Abs(trackat0.GetC()));
4059 if ( bestLong->GetNumberOfClusters()>3 ){
4060 dist[itsindex] = trackat0.GetY();
4061 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4062 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4063 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4064 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4066 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
4067 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
4068 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
4070 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
4071 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
4075 if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
4076 dist[itsindex] = bestConst->GetD(0);
4077 norm[itsindex] = bestConst->GetDnorm(0);
4078 normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4079 normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4080 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4082 dist[itsindex] = trackat0.GetY();
4083 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4084 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4085 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4086 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4087 if (TMath::Abs(trackat0.GetTgl())>1.05){
4088 if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
4089 if (normdist[itsindex]>3) {
4090 minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
4096 //-----------------------------------------------------------
4097 //Forbid primary track candidates -
4099 //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
4100 //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
4101 //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
4102 //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
4103 //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
4104 //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
4105 //-----------------------------------------------------------
4107 if (bestLong->GetNumberOfClusters()<4 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5) forbidden[itsindex]=kTRUE;
4108 if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5) forbidden[itsindex]=kTRUE;
4109 if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>=0 && bestConst->GetClIndex(1)>=0 ) forbidden[itsindex]=kTRUE;
4110 if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>=0) forbidden[itsindex]=kTRUE;
4111 if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2) forbidden[itsindex]=kTRUE;
4112 if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1) forbidden[itsindex]=kTRUE;
4113 if (bestConst->GetNormChi2(0)<2.5) {
4114 minPointAngle[itsindex]= 0.9999;
4115 maxr[itsindex] = 10;
4119 //forbid daughter kink candidates
4121 if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
4122 Bool_t isElectron = kTRUE;
4123 Bool_t isProton = kTRUE;
4125 esdtrack->GetESDpid(pid);
4126 for (Int_t i=1;i<5;i++){
4127 if (pid[0]<pid[i]) isElectron= kFALSE;
4128 if (pid[4]<pid[i]) isProton= kFALSE;
4131 forbidden[itsindex]=kFALSE;
4132 normdist[itsindex]*=-1;
4135 if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
4136 normdist[itsindex]*=-1;
4140 // Causality cuts in TPC volume
4142 if (esdtrack->GetTPCdensity(0,10) >0.6) maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
4143 if (esdtrack->GetTPCdensity(10,30)>0.6) maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
4144 if (esdtrack->GetTPCdensity(20,40)>0.6) maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
4145 if (esdtrack->GetTPCdensity(30,50)>0.6) maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
4147 if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
4153 "Tr1.="<<((bestConst)? bestConst:dummy)<<
4155 "Tr3.="<<&trackat0<<
4157 "Dist="<<dist[itsindex]<<
4158 "ND0="<<normdist0[itsindex]<<
4159 "ND1="<<normdist1[itsindex]<<
4160 "ND="<<normdist[itsindex]<<
4161 "Pz="<<primvertex[2]<<
4162 "Forbid="<<forbidden[itsindex]<<
4166 trackarray.AddAt(best,itsindex);
4167 trackarrayc.AddAt(bestConst,itsindex);
4168 trackarrayl.AddAt(bestLong,itsindex);
4169 new (&helixes[itsindex]) AliHelix(*best);
4174 // first iterration of V0 finder
4176 for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
4177 Int_t itrack0 = itsmap[iesd0];
4178 if (forbidden[itrack0]) continue;
4179 AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
4180 if (!btrack0) continue;
4181 if (btrack0->GetSign()>0 && !AliITSReconstructor::GetRecoParam()->GetStoreLikeSignV0s()) continue;
4182 AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
4184 for (Int_t iesd1=iesd0+1;iesd1<ntracks;iesd1++){
4185 Int_t itrack1 = itsmap[iesd1];
4186 if (forbidden[itrack1]) continue;
4188 AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
4189 if (!btrack1) continue;
4190 if (btrack1->GetSign()<0 && !AliITSReconstructor::GetRecoParam()->GetStoreLikeSignV0s()) continue;
4191 Bool_t isGold = kFALSE;
4192 if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
4195 AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
4196 AliHelix &h1 = helixes[itrack0];
4197 AliHelix &h2 = helixes[itrack1];
4199 // find linear distance
4204 Double_t phase[2][2],radius[2];
4205 Int_t points = h1.GetRPHIintersections(h2, phase, radius);
4206 if (points==0) continue;
4207 Double_t delta[2]={1000000,1000000};
4209 h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
4211 if (radius[1]<rmin) rmin = radius[1];
4212 h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
4214 rmin = TMath::Sqrt(rmin);
4215 Double_t distance = 0;
4216 Double_t radiusC = 0;
4218 if (points==1 || delta[0]<delta[1]){
4219 distance = TMath::Sqrt(delta[0]);
4220 radiusC = TMath::Sqrt(radius[0]);
4222 distance = TMath::Sqrt(delta[1]);
4223 radiusC = TMath::Sqrt(radius[1]);
4226 if (radiusC<TMath::Max(minr[itrack0],minr[itrack1])) continue;
4227 if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1])) continue;
4228 Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
4229 if (distance>maxDist) continue;
4230 Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
4231 if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
4234 // Double_t distance = TestV0(h1,h2,pvertex,rmin);
4236 // if (distance>maxDist) continue;
4237 // if (pvertex->GetRr()<kMinR) continue;
4238 // if (pvertex->GetRr()>kMaxR) continue;
4239 AliITStrackMI * track0=btrack0;
4240 AliITStrackMI * track1=btrack1;
4241 // if (pvertex->GetRr()<3.5){
4243 //use longest tracks inside the pipe
4244 track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4245 track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4249 pvertex->SetParamN(*track0);
4250 pvertex->SetParamP(*track1);
4251 pvertex->Update(primvertex);
4252 pvertex->SetClusters(track0->ClIndex(),track1->ClIndex()); // register clusters
4254 if (pvertex->GetRr()<kMinR) continue;
4255 if (pvertex->GetRr()>kMaxR) continue;
4256 if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4257 //Bo: if (pvertex->GetDist2()>maxDist) continue;
4258 if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4259 //Bo: pvertex->SetLab(0,track0->GetLabel());
4260 //Bo: pvertex->SetLab(1,track1->GetLabel());
4261 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4262 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4264 AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4265 AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4269 TObjArray * array0b = (TObjArray*)fBestHypothesys.At(itrack0);
4270 if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4271 fCurrentEsdTrack = itrack0;
4272 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4274 TObjArray * array1b = (TObjArray*)fBestHypothesys.At(itrack1);
4275 if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4276 fCurrentEsdTrack = itrack1;
4277 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4280 AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4281 AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4282 AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4283 AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4285 Float_t minchi2before0=16;
4286 Float_t minchi2before1=16;
4287 Float_t minchi2after0 =16;
4288 Float_t minchi2after1 =16;
4289 Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
4290 Int_t maxLayer = GetNearestLayer(xrp); //I.B.
4292 if (array0b) for (Int_t i=0;i<5;i++){
4293 // best track after vertex
4294 AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4295 if (!btrack) continue;
4296 if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4297 // if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4298 if (btrack->GetX()<pvertex->GetRr()-2.) {
4299 if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4302 if (maxLayer<3){ // take prim vertex as additional measurement
4303 if (normdist[itrack0]>0 && htrackc0){
4304 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4306 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4310 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4312 if (btrack->GetClIndex(ilayer)<0){
4316 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4317 for (Int_t itrack=0;itrack<4;itrack++){
4318 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4319 sumchi2+=18.; //shared cluster
4323 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4324 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4328 if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4330 continue; //safety space - Geo manager will give exact layer
4333 minchi2after0 = btrack->GetNormChi2(i);
4336 if (array1b) for (Int_t i=0;i<5;i++){
4337 // best track after vertex
4338 AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4339 if (!btrack) continue;
4340 if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4341 // if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4342 if (btrack->GetX()<pvertex->GetRr()-2){
4343 if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4346 if (maxLayer<3){ // take prim vertex as additional measurement
4347 if (normdist[itrack1]>0 && htrackc1){
4348 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4350 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4354 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4356 if (btrack->GetClIndex(ilayer)<0){
4360 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4361 for (Int_t itrack=0;itrack<4;itrack++){
4362 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4363 sumchi2+=18.; //shared cluster
4367 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4368 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4372 if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4374 continue; //safety space - Geo manager will give exact layer
4377 minchi2after1 = btrack->GetNormChi2(i);
4381 // position resolution - used for DCA cut
4382 Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4383 (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4384 (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4385 sigmad =TMath::Sqrt(sigmad)+0.04;
4386 if (pvertex->GetRr()>50){
4387 Double_t cov0[15],cov1[15];
4388 track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4389 track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4390 sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4391 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4392 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4393 sigmad =TMath::Sqrt(sigmad)+0.05;
4397 vertex2.SetParamN(*track0b);
4398 vertex2.SetParamP(*track1b);
4399 vertex2.Update(primvertex);
4400 //Bo: if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4401 if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4402 pvertex->SetParamN(*track0b);
4403 pvertex->SetParamP(*track1b);
4404 pvertex->Update(primvertex);
4405 pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex()); // register clusters
4406 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4407 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4409 pvertex->SetDistSigma(sigmad);
4410 //Bo: pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4411 pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4413 // define likelihhod and causalities
4415 Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4417 Float_t fnorm0 = normdist[itrack0];
4418 if (fnorm0<0) fnorm0*=-3;
4419 Float_t fnorm1 = normdist[itrack1];
4420 if (fnorm1<0) fnorm1*=-3;
4421 if ((pvertex->GetAnglep()[2]>0.1) || ( (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 ) || (pvertex->GetRr()<3)){
4422 pb0 = TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4423 pb1 = TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4425 pvertex->SetChi2Before(normdist[itrack0]);
4426 pvertex->SetChi2After(normdist[itrack1]);
4427 pvertex->SetNAfter(0);
4428 pvertex->SetNBefore(0);
4430 pvertex->SetChi2Before(minchi2before0);
4431 pvertex->SetChi2After(minchi2before1);
4432 if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4433 pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4434 pb1 = TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4436 pvertex->SetNAfter(maxLayer);
4437 pvertex->SetNBefore(maxLayer);
4439 if (pvertex->GetRr()<90){
4440 pa0 *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4441 pa1 *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4443 if (pvertex->GetRr()<20){
4444 pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4445 pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4448 pvertex->SetCausality(pb0,pb1,pa0,pa1);
4450 // Likelihood calculations - apply cuts
4452 Bool_t v0OK = kTRUE;
4453 Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4454 p12 += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4455 p12 = TMath::Sqrt(p12); // "mean" momenta
4456 Float_t sigmap0 = 0.0001+0.001/(0.1+pvertex->GetRr());
4457 Float_t sigmap = 0.5*sigmap0*(0.6+0.4*p12); // "resolution: of point angle - as a function of radius and momenta
4459 Float_t causalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4460 Float_t causalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4461 TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4463 //Bo: Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4464 Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4465 Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4467 Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4468 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4469 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4470 0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4472 if (causalityA<kCausality0Cut) v0OK = kFALSE;
4473 if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut) v0OK = kFALSE;
4474 if (likelihood1<kLikelihood1Cut) v0OK = kFALSE;
4475 if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4480 Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4482 "Tr0.="<<track0<< //best without constrain
4483 "Tr1.="<<track1<< //best without constrain
4484 "Tr0B.="<<track0b<< //best without constrain after vertex
4485 "Tr1B.="<<track1b<< //best without constrain after vertex
4486 "Tr0C.="<<htrackc0<< //best with constrain if exist
4487 "Tr1C.="<<htrackc1<< //best with constrain if exist
4488 "Tr0L.="<<track0l<< //longest best
4489 "Tr1L.="<<track1l<< //longest best
4490 "Esd0.="<<track0->GetESDtrack()<< // esd track0 params
4491 "Esd1.="<<track1->GetESDtrack()<< // esd track1 params
4492 "V0.="<<pvertex<< //vertex properties
4493 "V0b.="<<&vertex2<< //vertex properties at "best" track
4494 "ND0="<<normdist[itrack0]<< //normalize distance for track0
4495 "ND1="<<normdist[itrack1]<< //normalize distance for track1
4497 // "RejectBase="<<rejectBase<< //rejection in First itteration
4503 //if (rejectBase) continue;
4505 pvertex->SetStatus(0);
4506 // if (rejectBase) {
4507 // pvertex->SetStatus(-100);
4509 if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4510 //Bo: pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4511 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4512 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4514 // AliV0vertex vertexjuri(*track0,*track1);
4515 // vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4516 // event->AddV0(&vertexjuri);
4517 pvertex->SetStatus(100);
4519 pvertex->SetOnFlyStatus(kTRUE);
4520 pvertex->ChangeMassHypothesis(kK0Short);
4521 event->AddV0(pvertex);
4527 // delete temporary arrays
4530 delete[] minPointAngle;
4542 //------------------------------------------------------------------------
4543 void AliITStrackerMI::RefitV02(const AliESDEvent *event)
4546 //try to refit V0s in the third path of the reconstruction
4548 TTreeSRedirector &cstream = *fDebugStreamer;
4550 Int_t nv0s = event->GetNumberOfV0s();
4551 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4553 for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4554 AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4555 if (!v0mi) continue;
4556 Int_t itrack0 = v0mi->GetIndex(0);
4557 Int_t itrack1 = v0mi->GetIndex(1);
4558 AliESDtrack *esd0 = event->GetTrack(itrack0);
4559 AliESDtrack *esd1 = event->GetTrack(itrack1);
4560 if (!esd0||!esd1) continue;
4561 AliITStrackMI tpc0(*esd0);
4562 AliITStrackMI tpc1(*esd1);
4563 Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4564 Double_t alpha =TMath::ATan2(y,x); //I.B.
4565 if (v0mi->GetRr()>85){
4566 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4567 v0temp.SetParamN(tpc0);
4568 v0temp.SetParamP(tpc1);
4569 v0temp.Update(primvertex);
4570 if (kFALSE) cstream<<"Refit"<<
4572 "V0refit.="<<&v0temp<<
4576 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4577 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4578 v0mi->SetParamN(tpc0);
4579 v0mi->SetParamP(tpc1);
4580 v0mi->Update(primvertex);
4585 if (v0mi->GetRr()>35){
4586 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4587 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4588 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4589 v0temp.SetParamN(tpc0);
4590 v0temp.SetParamP(tpc1);
4591 v0temp.Update(primvertex);
4592 if (kFALSE) cstream<<"Refit"<<
4594 "V0refit.="<<&v0temp<<
4598 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4599 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4600 v0mi->SetParamN(tpc0);
4601 v0mi->SetParamP(tpc1);
4602 v0mi->Update(primvertex);
4607 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4608 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4609 // if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4610 if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4611 v0temp.SetParamN(tpc0);
4612 v0temp.SetParamP(tpc1);
4613 v0temp.Update(primvertex);
4614 if (kFALSE) cstream<<"Refit"<<
4616 "V0refit.="<<&v0temp<<
4620 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4621 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4622 v0mi->SetParamN(tpc0);
4623 v0mi->SetParamP(tpc1);
4624 v0mi->Update(primvertex);
4630 //------------------------------------------------------------------------
4631 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4632 //--------------------------------------------------------------------
4633 // Fill a look-up table with mean material
4634 //--------------------------------------------------------------------
4638 Double_t point1[3],point2[3];
4639 Double_t phi,cosphi,sinphi,z;
4640 // 0-5 layers, 6 pipe, 7-8 shields
4641 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4642 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4644 Int_t ifirst=0,ilast=0;
4645 if(material.Contains("Pipe")) {
4647 } else if(material.Contains("Shields")) {
4649 } else if(material.Contains("Layers")) {
4652 Error("BuildMaterialLUT","Wrong layer name\n");
4655 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4656 Double_t param[5]={0.,0.,0.,0.,0.};
4657 for (Int_t i=0; i<n; i++) {
4658 phi = 2.*TMath::Pi()*gRandom->Rndm();
4659 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4660 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4661 point1[0] = rmin[imat]*cosphi;
4662 point1[1] = rmin[imat]*sinphi;
4664 point2[0] = rmax[imat]*cosphi;
4665 point2[1] = rmax[imat]*sinphi;
4667 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4668 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4670 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4672 fxOverX0Layer[imat] = param[1];
4673 fxTimesRhoLayer[imat] = param[0]*param[4];
4674 } else if(imat==6) {
4675 fxOverX0Pipe = param[1];
4676 fxTimesRhoPipe = param[0]*param[4];
4677 } else if(imat==7) {
4678 fxOverX0Shield[0] = param[1];
4679 fxTimesRhoShield[0] = param[0]*param[4];
4680 } else if(imat==8) {
4681 fxOverX0Shield[1] = param[1];
4682 fxTimesRhoShield[1] = param[0]*param[4];
4686 printf("%s\n",material.Data());
4687 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4688 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4689 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4690 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4691 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4692 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4693 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4694 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4695 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4699 //------------------------------------------------------------------------
4700 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4701 TString direction) {
4702 //-------------------------------------------------------------------
4703 // Propagate beyond beam pipe and correct for material
4704 // (material budget in different ways according to fUseTGeo value)
4705 // Add time if going outward (PropagateTo or PropagateToTGeo)
4706 //-------------------------------------------------------------------
4708 // Define budget mode:
4709 // 0: material from AliITSRecoParam (hard coded)
4710 // 1: material from TGeo in one step (on the fly)
4711 // 2: material from lut
4712 // 3: material from TGeo in one step (same for all hypotheses)
4725 if(fTrackingPhase.Contains("Clusters2Tracks"))
4726 { mode=3; } else { mode=1; }
4729 if(fTrackingPhase.Contains("Clusters2Tracks"))
4730 { mode=3; } else { mode=2; }
4736 if(fTrackingPhase.Contains("Default")) mode=0;
4738 Int_t index=fCurrentEsdTrack;
4740 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4741 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4743 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4745 Double_t xOverX0,x0,lengthTimesMeanDensity;
4749 xOverX0 = AliITSRecoParam::GetdPipe();
4750 x0 = AliITSRecoParam::GetX0Be();
4751 lengthTimesMeanDensity = xOverX0*x0;
4752 lengthTimesMeanDensity *= dir;
4753 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4756 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4759 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4760 xOverX0 = fxOverX0Pipe;
4761 lengthTimesMeanDensity = fxTimesRhoPipe;
4762 lengthTimesMeanDensity *= dir;
4763 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4766 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4767 if(fxOverX0PipeTrks[index]<0) {
4768 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4769 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4770 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4771 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4772 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4775 xOverX0 = fxOverX0PipeTrks[index];
4776 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4777 lengthTimesMeanDensity *= dir;
4778 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4784 //------------------------------------------------------------------------
4785 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4787 TString direction) {
4788 //-------------------------------------------------------------------
4789 // Propagate beyond SPD or SDD shield and correct for material
4790 // (material budget in different ways according to fUseTGeo value)
4791 // Add time if going outward (PropagateTo or PropagateToTGeo)
4792 //-------------------------------------------------------------------
4794 // Define budget mode:
4795 // 0: material from AliITSRecoParam (hard coded)
4796 // 1: material from TGeo in steps of X cm (on the fly)
4797 // X = AliITSRecoParam::GetStepSizeTGeo()
4798 // 2: material from lut
4799 // 3: material from TGeo in one step (same for all hypotheses)
4812 if(fTrackingPhase.Contains("Clusters2Tracks"))
4813 { mode=3; } else { mode=1; }
4816 if(fTrackingPhase.Contains("Clusters2Tracks"))
4817 { mode=3; } else { mode=2; }
4823 if(fTrackingPhase.Contains("Default")) mode=0;
4825 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4827 Int_t shieldindex=0;
4828 if (shield.Contains("SDD")) { // SDDouter
4829 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4831 } else if (shield.Contains("SPD")) { // SPDouter
4832 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4835 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4839 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4841 Int_t index=2*fCurrentEsdTrack+shieldindex;
4843 Double_t xOverX0,x0,lengthTimesMeanDensity;
4848 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4849 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4850 lengthTimesMeanDensity = xOverX0*x0;
4851 lengthTimesMeanDensity *= dir;
4852 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4855 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4856 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4859 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4860 xOverX0 = fxOverX0Shield[shieldindex];
4861 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4862 lengthTimesMeanDensity *= dir;
4863 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4866 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4867 if(fxOverX0ShieldTrks[index]<0) {
4868 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4869 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4870 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4871 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4872 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4875 xOverX0 = fxOverX0ShieldTrks[index];
4876 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4877 lengthTimesMeanDensity *= dir;
4878 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4884 //------------------------------------------------------------------------
4885 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4887 Double_t oldGlobXYZ[3],
4888 TString direction) {
4889 //-------------------------------------------------------------------
4890 // Propagate beyond layer and correct for material
4891 // (material budget in different ways according to fUseTGeo value)
4892 // Add time if going outward (PropagateTo or PropagateToTGeo)
4893 //-------------------------------------------------------------------
4895 // Define budget mode:
4896 // 0: material from AliITSRecoParam (hard coded)
4897 // 1: material from TGeo in stepsof X cm (on the fly)
4898 // X = AliITSRecoParam::GetStepSizeTGeo()
4899 // 2: material from lut
4900 // 3: material from TGeo in one step (same for all hypotheses)
4913 if(fTrackingPhase.Contains("Clusters2Tracks"))
4914 { mode=3; } else { mode=1; }
4917 if(fTrackingPhase.Contains("Clusters2Tracks"))
4918 { mode=3; } else { mode=2; }
4924 if(fTrackingPhase.Contains("Default")) mode=0;
4926 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4928 Double_t r=fgLayers[layerindex].GetR();
4929 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4931 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4933 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4935 Int_t index=6*fCurrentEsdTrack+layerindex;
4938 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4941 // back before material (no correction)
4943 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4944 if (!t->GetLocalXat(rOld,xOld)) return 0;
4945 if (!t->Propagate(xOld)) return 0;
4949 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4950 lengthTimesMeanDensity = xOverX0*x0;
4951 lengthTimesMeanDensity *= dir;
4952 // Bring the track beyond the material
4953 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4956 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4957 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4960 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4961 xOverX0 = fxOverX0Layer[layerindex];
4962 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4963 lengthTimesMeanDensity *= dir;
4964 // Bring the track beyond the material
4965 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4968 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4969 if(fxOverX0LayerTrks[index]<0) {
4970 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4971 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4972 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4973 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4974 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4975 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4978 xOverX0 = fxOverX0LayerTrks[index];
4979 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4980 lengthTimesMeanDensity *= dir;
4981 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4988 //------------------------------------------------------------------------
4989 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4990 //-----------------------------------------------------------------
4991 // Initialize LUT for storing material for each prolonged track
4992 //-----------------------------------------------------------------
4993 fxOverX0PipeTrks = new Float_t[ntracks];
4994 fxTimesRhoPipeTrks = new Float_t[ntracks];
4995 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4996 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4997 fxOverX0LayerTrks = new Float_t[ntracks*6];
4998 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
5000 for(Int_t i=0; i<ntracks; i++) {
5001 fxOverX0PipeTrks[i] = -1.;
5002 fxTimesRhoPipeTrks[i] = -1.;
5004 for(Int_t j=0; j<ntracks*2; j++) {
5005 fxOverX0ShieldTrks[j] = -1.;
5006 fxTimesRhoShieldTrks[j] = -1.;
5008 for(Int_t k=0; k<ntracks*6; k++) {
5009 fxOverX0LayerTrks[k] = -1.;
5010 fxTimesRhoLayerTrks[k] = -1.;
5017 //------------------------------------------------------------------------
5018 void AliITStrackerMI::DeleteTrksMaterialLUT() {
5019 //-----------------------------------------------------------------
5020 // Delete LUT for storing material for each prolonged track
5021 //-----------------------------------------------------------------
5022 if(fxOverX0PipeTrks) {
5023 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
5025 if(fxOverX0ShieldTrks) {
5026 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
5029 if(fxOverX0LayerTrks) {
5030 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
5032 if(fxTimesRhoPipeTrks) {
5033 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
5035 if(fxTimesRhoShieldTrks) {
5036 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
5038 if(fxTimesRhoLayerTrks) {
5039 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
5043 //------------------------------------------------------------------------
5044 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
5045 Int_t ilayer,Int_t idet) const {
5046 //-----------------------------------------------------------------
5047 // This method is used to decide whether to allow a prolongation
5048 // without clusters, because we want to skip the layer.
5049 // In this case the return value is > 0:
5050 // return 1: the user requested to skip a layer
5051 // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
5052 //-----------------------------------------------------------------
5054 if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
5056 if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
5057 // check if track will cross SPD outer layer
5058 Double_t phiAtSPD2,zAtSPD2;
5059 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
5060 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
5066 //------------------------------------------------------------------------
5067 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
5068 Int_t ilayer,Int_t idet,
5069 Double_t dz,Double_t dy,
5070 Bool_t noClusters) const {
5071 //-----------------------------------------------------------------
5072 // This method is used to decide whether to allow a prolongation
5073 // without clusters, because there is a dead zone in the road.
5074 // In this case the return value is > 0:
5075 // return 1: dead zone at z=0,+-7cm in SPD
5076 // return 2: all road is "bad" (dead or noisy) from the OCDB
5077 // return 3: something "bad" (dead or noisy) from the OCDB
5078 //-----------------------------------------------------------------
5080 // check dead zones at z=0,+-7cm in the SPD
5081 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
5082 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5083 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5084 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
5085 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5086 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5087 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
5088 for (Int_t i=0; i<3; i++)
5089 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
5090 AliDebug(2,Form("crack SPD %d",ilayer));
5095 // check bad zones from OCDB
5096 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
5098 if (idet<0) return 0;
5100 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
5103 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
5104 if (ilayer==0 || ilayer==1) { // ---------- SPD
5106 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
5108 detSizeFactorX *= 2.;
5109 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
5112 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
5113 if (detType==2) segm->SetLayer(ilayer+1);
5114 Float_t detSizeX = detSizeFactorX*segm->Dx();
5115 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
5117 // check if the road overlaps with bad chips
5119 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
5120 Float_t zlocmin = zloc-dz;
5121 Float_t zlocmax = zloc+dz;
5122 Float_t xlocmin = xloc-dy;
5123 Float_t xlocmax = xloc+dy;
5124 Int_t chipsInRoad[100];
5126 // check if road goes out of detector
5127 Bool_t touchNeighbourDet=kFALSE;
5128 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
5129 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
5130 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5131 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5132 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));
5134 // check if this detector is bad
5136 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
5137 if(!touchNeighbourDet) {
5138 return 2; // all detectors in road are bad
5140 return 3; // at least one is bad
5144 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
5145 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
5146 if (!nChipsInRoad) return 0;
5148 Bool_t anyBad=kFALSE,anyGood=kFALSE;
5149 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
5150 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
5151 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
5152 if (det.IsChipBad(chipsInRoad[iCh])) {
5160 if(!touchNeighbourDet) {
5161 AliDebug(2,"all bad in road");
5162 return 2; // all chips in road are bad
5164 return 3; // at least a bad chip in road
5169 AliDebug(2,"at least a bad in road");
5170 return 3; // at least a bad chip in road
5174 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
5175 || ilayer==4 || ilayer==5 // SSD
5176 || !noClusters) return 0;
5178 // There are no clusters in road: check if there is at least
5179 // a bad SPD pixel or SDD anode
5181 Int_t idetInITS=idet;
5182 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
5184 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
5185 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
5188 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
5192 //------------------------------------------------------------------------
5193 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
5194 const AliITStrackMI *track,
5195 Float_t &xloc,Float_t &zloc) const {
5196 //-----------------------------------------------------------------
5197 // Gives position of track in local module ref. frame
5198 //-----------------------------------------------------------------
5203 if(idet<0) return kFALSE;
5205 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
5207 Int_t lad = Int_t(idet/ndet) + 1;
5209 Int_t det = idet - (lad-1)*ndet + 1;
5211 Double_t xyzGlob[3],xyzLoc[3];
5213 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
5214 // take into account the misalignment: xyz at real detector plane
5215 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
5217 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
5219 xloc = (Float_t)xyzLoc[0];
5220 zloc = (Float_t)xyzLoc[2];
5224 //------------------------------------------------------------------------
5225 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
5227 // Method to be optimized further:
5228 // Aim: decide whether a track can be used for PlaneEff evaluation
5229 // the decision is taken based on the track quality at the layer under study
5230 // no information on the clusters on this layer has to be used
5231 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
5232 // the cut is done on number of sigmas from the boundaries
5234 // Input: Actual track, layer [0,5] under study
5236 // Return: kTRUE if this is a good track
5238 // it will apply a pre-selection to obtain good quality tracks.
5239 // Here also you will have the possibility to put a control on the
5240 // impact point of the track on the basic block, in order to exclude border regions
5241 // this will be done by calling a proper method of the AliITSPlaneEff class.
5243 // input: AliITStrackMI* track, ilayer= layer number [0,5]
5244 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
5246 Int_t index[AliITSgeomTGeo::kNLayers];
5248 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
5250 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
5251 index[k]=clusters[k];
5255 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
5256 AliITSlayer &layer=fgLayers[ilayer];
5257 Double_t r=layer.GetR();
5258 AliITStrackMI tmp(*track);
5260 // require a minimal number of cluster in other layers and eventually clusters in closest layers
5262 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
5263 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
5264 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
5265 if (tmp.GetClIndex(lay)>=0) ncl++;
5267 Bool_t nextout = kFALSE;
5268 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
5269 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
5270 Bool_t nextin = kFALSE;
5271 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
5272 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
5273 if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
5275 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
5276 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
5277 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
5278 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
5282 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
5283 Int_t idet=layer.FindDetectorIndex(phi,z);
5284 if(idet<0) { AliInfo(Form("cannot find detector"));
5287 // here check if it has good Chi Square.
5289 //propagate to the intersection with the detector plane
5290 const AliITSdetector &det=layer.GetDetector(idet);
5291 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
5295 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
5296 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5297 if(key>fPlaneEff->Nblock()) return kFALSE;
5298 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
5299 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
5301 // DEFINITION OF SEARCH ROAD FOR accepting a track
5303 //For the time being they are hard-wired, later on from AliITSRecoParam
5304 // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
5305 // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
5308 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5309 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5310 // done for RecPoints
5312 // exclude tracks at boundary between detectors
5313 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
5314 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
5315 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
5316 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
5317 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5319 if ( (locx-dx < blockXmn+boundaryWidth) ||
5320 (locx+dx > blockXmx-boundaryWidth) ||
5321 (locz-dz < blockZmn+boundaryWidth) ||
5322 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5325 //------------------------------------------------------------------------
5326 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
5328 // This Method has to be optimized! For the time-being it uses the same criteria
5329 // as those used in the search of extra clusters for overlapping modules.
5331 // Method Purpose: estabilish whether a track has produced a recpoint or not
5332 // in the layer under study (For Plane efficiency)
5334 // inputs: AliITStrackMI* track (pointer to a usable track)
5336 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5337 // traversed by this very track. In details:
5338 // - if a cluster can be associated to the track then call
5339 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5340 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5343 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5344 AliITSlayer &layer=fgLayers[ilayer];
5345 Double_t r=layer.GetR();
5346 AliITStrackMI tmp(*track);
5350 if (!tmp.GetPhiZat(r,phi,z)) return;
5351 Int_t idet=layer.FindDetectorIndex(phi,z);
5353 if(idet<0) { AliInfo(Form("cannot find detector"));
5357 //propagate to the intersection with the detector plane
5358 const AliITSdetector &det=layer.GetDetector(idet);
5359 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5363 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5365 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5366 TMath::Sqrt(tmp.GetSigmaZ2() +
5367 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5368 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5369 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5370 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5371 TMath::Sqrt(tmp.GetSigmaY2() +
5372 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5373 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5374 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5376 // road in global (rphi,z) [i.e. in tracking ref. system]
5377 Double_t zmin = tmp.GetZ() - dz;
5378 Double_t zmax = tmp.GetZ() + dz;
5379 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5380 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5382 // select clusters in road
5383 layer.SelectClusters(zmin,zmax,ymin,ymax);
5385 // Define criteria for track-cluster association
5386 Double_t msz = tmp.GetSigmaZ2() +
5387 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5388 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5389 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5390 Double_t msy = tmp.GetSigmaY2() +
5391 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5392 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5393 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5394 if (tmp.GetConstrain()) {
5395 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5396 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5398 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5399 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5401 msz = 1./msz; // 1/RoadZ^2
5402 msy = 1./msy; // 1/RoadY^2
5405 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5407 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5408 //Double_t tolerance=0.2;
5409 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5410 idetc = cl->GetDetectorIndex();
5411 if(idet!=idetc) continue;
5412 //Int_t ilay = cl->GetLayer();
5414 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5415 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5417 Double_t chi2=tmp.GetPredictedChi2(cl);
5418 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5422 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5424 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5425 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5426 if(key>fPlaneEff->Nblock()) return;
5427 Bool_t found=kFALSE;
5430 while ((cl=layer.GetNextCluster(clidx))!=0) {
5431 idetc = cl->GetDetectorIndex();
5432 if(idet!=idetc) continue;
5433 // here real control to see whether the cluster can be associated to the track.
5434 // cluster not associated to track
5435 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5436 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5437 // calculate track-clusters chi2
5438 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5439 // in particular, the error associated to the cluster
5440 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5442 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5444 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5445 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5446 // track->SetExtraModule(ilayer,idetExtra);
5448 if(!fPlaneEff->UpDatePlaneEff(found,key))
5449 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5450 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5451 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5452 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5453 Int_t cltype[2]={-999,-999};
5457 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5458 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5461 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5462 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5463 cltype[0]=layer.GetCluster(ci)->GetNy();
5464 cltype[1]=layer.GetCluster(ci)->GetNz();
5466 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5467 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5468 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5469 // It is computed properly by calling the method
5470 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5472 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5473 //if (tmp.PropagateTo(x,0.,0.)) {
5474 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5475 AliCluster c(*layer.GetCluster(ci));
5476 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5477 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5478 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5479 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5480 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5483 fPlaneEff->FillHistos(key,found,tr,clu,cltype);