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 // dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23 // Params moved to AliITSRecoParam by: Andrea Dainese, INFN
24 // Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
25 //-------------------------------------------------------------------------
29 #include <TTreeStream.h>
30 #include <TDatabasePDG.h>
35 #include "AliESDEvent.h"
36 #include "AliESDtrack.h"
37 #include "AliESDVertex.h"
40 #include "AliITSRecPoint.h"
41 #include "AliITSgeomTGeo.h"
42 #include "AliITSReconstructor.h"
43 #include "AliTrackPointArray.h"
44 #include "AliAlignObj.h"
45 #include "AliITSClusterParam.h"
46 #include "AliITSPlaneEff.h"
47 #include "AliCDBManager.h"
48 #include "AliCDBEntry.h"
49 #include "AliITSCalibrationSPD.h"
50 #include "AliITSCalibrationSDD.h"
51 #include "AliITSCalibrationSSD.h"
52 #include "AliITStrackerMI.h"
54 ClassImp(AliITStrackerMI)
56 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
58 AliITStrackerMI::AliITStrackerMI():AliTracker(),
68 fLastLayerToTrackTo(0),
71 fTrackingPhase("Default"),
77 fxTimesRhoPipeTrks(0),
78 fxOverX0ShieldTrks(0),
79 fxTimesRhoShieldTrks(0),
81 fxTimesRhoLayerTrks(0),
86 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
87 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
88 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
90 //------------------------------------------------------------------------
91 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
92 fI(AliITSgeomTGeo::GetNLayers()),
101 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
104 fTrackingPhase("Default"),
110 fxTimesRhoPipeTrks(0),
111 fxOverX0ShieldTrks(0),
112 fxTimesRhoShieldTrks(0),
113 fxOverX0LayerTrks(0),
114 fxTimesRhoLayerTrks(0),
117 //--------------------------------------------------------------------
118 //This is the AliITStrackerMI constructor
119 //--------------------------------------------------------------------
121 AliWarning("\"geom\" is actually a dummy argument !");
127 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
128 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
129 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
131 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
132 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
133 Double_t poff=TMath::ATan2(y,x);
135 Double_t r=TMath::Sqrt(x*x + y*y);
137 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
138 r += TMath::Sqrt(x*x + y*y);
139 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
140 r += TMath::Sqrt(x*x + y*y);
141 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
142 r += TMath::Sqrt(x*x + y*y);
145 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
147 for (Int_t j=1; j<nlad+1; j++) {
148 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
149 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
150 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
152 Double_t txyz[3]={0.}, xyz[3]={0.};
153 m.LocalToMaster(txyz,xyz);
154 Double_t r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
155 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
157 if (phi<0) phi+=TMath::TwoPi();
158 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
160 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
161 new(&det) AliITSdetector(r,phi);
167 fI=AliITSgeomTGeo::GetNLayers();
170 fConstraint[0]=1; fConstraint[1]=0;
172 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
173 AliITSReconstructor::GetRecoParam()->GetYVdef(),
174 AliITSReconstructor::GetRecoParam()->GetZVdef()};
175 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
176 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
177 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
178 SetVertex(xyzVtx,ersVtx);
180 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=AliITSRecoParam::GetLayersNotToSkip(i);
181 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
182 for (Int_t i=0;i<100000;i++){
183 fBestTrackIndex[i]=0;
186 // store positions of centre of SPD modules (in z)
188 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
189 fSPDdetzcentre[0] = tr[2];
190 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
191 fSPDdetzcentre[1] = tr[2];
192 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
193 fSPDdetzcentre[2] = tr[2];
194 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
195 fSPDdetzcentre[3] = tr[2];
197 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
198 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
199 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
203 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
204 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
206 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
209 //------------------------------------------------------------------------
210 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
212 fBestTrack(tracker.fBestTrack),
213 fTrackToFollow(tracker.fTrackToFollow),
214 fTrackHypothesys(tracker.fTrackHypothesys),
215 fBestHypothesys(tracker.fBestHypothesys),
216 fOriginal(tracker.fOriginal),
217 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
218 fPass(tracker.fPass),
219 fAfterV0(tracker.fAfterV0),
220 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
221 fCoefficients(tracker.fCoefficients),
223 fTrackingPhase(tracker.fTrackingPhase),
224 fUseTGeo(tracker.fUseTGeo),
225 fNtracks(tracker.fNtracks),
226 fxOverX0Pipe(tracker.fxOverX0Pipe),
227 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
229 fxTimesRhoPipeTrks(0),
230 fxOverX0ShieldTrks(0),
231 fxTimesRhoShieldTrks(0),
232 fxOverX0LayerTrks(0),
233 fxTimesRhoLayerTrks(0),
234 fDebugStreamer(tracker.fDebugStreamer),
235 fPlaneEff(tracker.fPlaneEff) {
239 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
242 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
243 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
246 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
247 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
250 //------------------------------------------------------------------------
251 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
252 //Assignment operator
253 this->~AliITStrackerMI();
254 new(this) AliITStrackerMI(tracker);
257 //------------------------------------------------------------------------
258 AliITStrackerMI::~AliITStrackerMI()
263 if (fCoefficients) delete [] fCoefficients;
264 DeleteTrksMaterialLUT();
265 if (fDebugStreamer) {
266 //fDebugStreamer->Close();
267 delete fDebugStreamer;
270 //------------------------------------------------------------------------
271 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
272 //--------------------------------------------------------------------
273 //This function set masks of the layers which must be not skipped
274 //--------------------------------------------------------------------
275 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=l[i];
277 //------------------------------------------------------------------------
278 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
279 //--------------------------------------------------------------------
280 //This function loads ITS clusters
281 //--------------------------------------------------------------------
282 TBranch *branch=cTree->GetBranch("ITSRecPoints");
284 Error("LoadClusters"," can't get the branch !\n");
288 TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
289 branch->SetAddress(&clusters);
293 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
294 Int_t ndet=fgLayers[i].GetNdetectors();
295 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
296 for (; j<jmax; j++) {
297 if (!cTree->GetEvent(j)) continue;
298 Int_t ncl=clusters->GetEntriesFast();
299 SignDeltas(clusters,GetZ());
302 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
303 detector=c->GetDetectorIndex();
305 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
307 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
310 // add dead zone "virtual" cluster in SPD, if there is a cluster within
311 // zwindow cm from the dead zone
312 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
313 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
314 Int_t lab[4] = {0,0,0,detector};
315 Int_t info[3] = {0,0,i};
316 Float_t q = 0.; // this identifies virtual clusters
317 Float_t hit[5] = {xdead,
319 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
320 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
322 Bool_t local = kTRUE;
323 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
324 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
325 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
326 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
327 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
328 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
329 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
330 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
331 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
332 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
333 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
334 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
335 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
336 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
337 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
338 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
339 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
340 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
341 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
343 } // "virtual" clusters in SPD
347 fgLayers[i].ResetRoad(); //road defined by the cluster density
348 fgLayers[i].SortClusters();
353 //------------------------------------------------------------------------
354 void AliITStrackerMI::UnloadClusters() {
355 //--------------------------------------------------------------------
356 //This function unloads ITS clusters
357 //--------------------------------------------------------------------
358 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
360 //------------------------------------------------------------------------
361 static Int_t CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
362 //--------------------------------------------------------------------
363 // Correction for the material between the TPC and the ITS
364 //--------------------------------------------------------------------
365 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
366 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
367 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
368 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
369 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
370 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
371 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
372 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
374 Error("CorrectForTPCtoITSDeadZoneMaterial","Track is already in the dead zone !");
380 //------------------------------------------------------------------------
381 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
382 //--------------------------------------------------------------------
383 // This functions reconstructs ITS tracks
384 // The clusters must be already loaded !
385 //--------------------------------------------------------------------
386 fTrackingPhase="Clusters2Tracks";
388 TObjArray itsTracks(15000);
390 fEsd = event; // store pointer to the esd
392 // temporary (for cosmics)
393 if(event->GetVertex()) {
394 TString title = event->GetVertex()->GetTitle();
395 if(title.Contains("cosmics")) {
396 Double_t xyz[3]={GetX(),GetY(),GetZ()};
397 Double_t exyz[3]={0.1,0.1,0.1};
403 {/* Read ESD tracks */
404 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
405 Int_t nentr=event->GetNumberOfTracks();
406 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
408 AliESDtrack *esd=event->GetTrack(nentr);
410 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
411 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
412 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
413 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
416 t=new AliITStrackMI(*esd);
417 } catch (const Char_t *msg) {
418 //Warning("Clusters2Tracks",msg);
422 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
423 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
426 // look at the ESD mass hypothesys !
427 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
429 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
431 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
432 //track - can be V0 according to TPC
434 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
438 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
442 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
446 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
451 t->SetReconstructed(kFALSE);
452 itsTracks.AddLast(t);
453 fOriginal.AddLast(t);
455 } /* End Read ESD tracks */
459 Int_t nentr=itsTracks.GetEntriesFast();
460 fTrackHypothesys.Expand(nentr);
461 fBestHypothesys.Expand(nentr);
462 MakeCoefficients(nentr);
463 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
465 // THE TWO TRACKING PASSES
466 for (fPass=0; fPass<2; fPass++) {
467 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
468 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
469 //cerr<<fPass<<" "<<fCurrentEsdTrack<<'\n';
470 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
471 if (t==0) continue; //this track has been already tracked
472 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
473 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
474 if (fConstraint[fPass]) {
475 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
476 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
479 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
481 ResetTrackToFollow(*t);
484 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
486 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
488 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
489 if (!besttrack) continue;
490 besttrack->SetLabel(tpcLabel);
491 // besttrack->CookdEdx();
493 besttrack->SetFakeRatio(1.);
494 CookLabel(besttrack,0.); //For comparison only
495 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
497 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
499 t->SetReconstructed(kTRUE);
502 GetBestHypothesysMIP(itsTracks);
503 } // end loop on the two tracking passes
505 //GetBestHypothesysMIP(itsTracks);
506 if(event->GetNumberOfV0s()>0) UpdateTPCV0(event);
507 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) FindV02(event);
509 //GetBestHypothesysMIP(itsTracks);
513 Int_t entries = fTrackHypothesys.GetEntriesFast();
514 for (Int_t ientry=0; ientry<entries; ientry++) {
515 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
516 if (array) array->Delete();
517 delete fTrackHypothesys.RemoveAt(ientry);
520 fTrackHypothesys.Delete();
521 fBestHypothesys.Delete();
523 delete [] fCoefficients;
525 DeleteTrksMaterialLUT();
527 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
529 fTrackingPhase="Default";
533 //------------------------------------------------------------------------
534 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
535 //--------------------------------------------------------------------
536 // This functions propagates reconstructed ITS tracks back
537 // The clusters must be loaded !
538 //--------------------------------------------------------------------
539 fTrackingPhase="PropagateBack";
540 Int_t nentr=event->GetNumberOfTracks();
541 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
544 for (Int_t i=0; i<nentr; i++) {
545 AliESDtrack *esd=event->GetTrack(i);
547 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
548 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
552 t=new AliITStrackMI(*esd);
553 } catch (const Char_t *msg) {
554 //Warning("PropagateBack",msg);
558 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
560 ResetTrackToFollow(*t);
562 // propagate to vertex [SR, GSI 17.02.2003]
563 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
564 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
565 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
566 fTrackToFollow.StartTimeIntegral();
567 // from vertex to outside pipe
568 CorrectForPipeMaterial(&fTrackToFollow,"outward");
571 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
572 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
573 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
577 fTrackToFollow.SetLabel(t->GetLabel());
578 //fTrackToFollow.CookdEdx();
579 CookLabel(&fTrackToFollow,0.); //For comparison only
580 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
581 //UseClusters(&fTrackToFollow);
587 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
589 fTrackingPhase="Default";
593 //------------------------------------------------------------------------
594 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
595 //--------------------------------------------------------------------
596 // This functions refits ITS tracks using the
597 // "inward propagated" TPC tracks
598 // The clusters must be loaded !
599 //--------------------------------------------------------------------
600 fTrackingPhase="RefitInward";
601 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) RefitV02(event);
602 Int_t nentr=event->GetNumberOfTracks();
603 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
606 for (Int_t i=0; i<nentr; i++) {
607 AliESDtrack *esd=event->GetTrack(i);
609 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
610 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
611 if (esd->GetStatus()&AliESDtrack::kTPCout)
612 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
616 t=new AliITStrackMI(*esd);
617 } catch (const Char_t *msg) {
618 //Warning("RefitInward",msg);
622 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
623 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
628 ResetTrackToFollow(*t);
629 fTrackToFollow.ResetClusters();
631 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
632 fTrackToFollow.ResetCovariance(10.);
635 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE)) {
636 fTrackToFollow.SetLabel(t->GetLabel());
637 // fTrackToFollow.CookdEdx();
638 CookdEdx(&fTrackToFollow);
640 CookLabel(&fTrackToFollow,0.0); //For comparison only
643 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
644 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
645 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
646 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
647 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
648 Float_t r[3]={0.,0.,0.};
650 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
657 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
659 fTrackingPhase="Default";
663 //------------------------------------------------------------------------
664 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
665 //--------------------------------------------------------------------
666 // Return pointer to a given cluster
667 //--------------------------------------------------------------------
668 Int_t l=(index & 0xf0000000) >> 28;
669 Int_t c=(index & 0x0fffffff) >> 00;
670 return fgLayers[l].GetCluster(c);
672 //------------------------------------------------------------------------
673 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
674 //--------------------------------------------------------------------
675 // Get track space point with index i
676 //--------------------------------------------------------------------
678 Int_t l=(index & 0xf0000000) >> 28;
679 Int_t c=(index & 0x0fffffff) >> 00;
680 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
681 Int_t idet = cl->GetDetectorIndex();
685 cl->GetGlobalXYZ(xyz);
686 cl->GetGlobalCov(cov);
689 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
692 iLayer = AliGeomManager::kSPD1;
695 iLayer = AliGeomManager::kSPD2;
698 iLayer = AliGeomManager::kSDD1;
701 iLayer = AliGeomManager::kSDD2;
704 iLayer = AliGeomManager::kSSD1;
707 iLayer = AliGeomManager::kSSD2;
710 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
713 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
714 p.SetVolumeID((UShort_t)volid);
717 //------------------------------------------------------------------------
718 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
719 AliTrackPoint& p, const AliESDtrack *t) {
720 //--------------------------------------------------------------------
721 // Get track space point with index i
722 // (assign error estimated during the tracking)
723 //--------------------------------------------------------------------
725 Int_t l=(index & 0xf0000000) >> 28;
726 Int_t c=(index & 0x0fffffff) >> 00;
727 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
728 Int_t idet = cl->GetDetectorIndex();
729 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
731 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
733 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
734 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
735 Double_t alpha = t->GetAlpha();
736 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
737 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,AliTracker::GetBz()));
738 phi += alpha-det.GetPhi();
739 Float_t tgphi = TMath::Tan(phi);
741 Float_t tgl = t->GetTgl(); // tgl about const along track
742 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
744 Float_t errlocalx,errlocalz;
745 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz);
749 cl->GetGlobalXYZ(xyz);
750 // cl->GetGlobalCov(cov);
751 Float_t pos[3] = {0.,0.,0.};
752 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
753 tmpcl.GetGlobalCov(cov);
757 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
760 iLayer = AliGeomManager::kSPD1;
763 iLayer = AliGeomManager::kSPD2;
766 iLayer = AliGeomManager::kSDD1;
769 iLayer = AliGeomManager::kSDD2;
772 iLayer = AliGeomManager::kSSD1;
775 iLayer = AliGeomManager::kSSD2;
778 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
781 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
782 p.SetVolumeID((UShort_t)volid);
785 //------------------------------------------------------------------------
786 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
788 //--------------------------------------------------------------------
789 // Follow prolongation tree
790 //--------------------------------------------------------------------
792 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
793 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
796 AliESDtrack * esd = otrack->GetESDtrack();
797 if (esd->GetV0Index(0)>0) {
798 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
799 // mapping of ESD track is different as ITS track in Containers
800 // Need something more stable
801 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
802 for (Int_t i=0;i<3;i++){
803 Int_t index = esd->GetV0Index(i);
805 AliESDv0 * vertex = fEsd->GetV0(index);
806 if (vertex->GetStatus()<0) continue; // rejected V0
808 if (esd->GetSign()>0) {
809 vertex->SetIndex(0,esdindex);
811 vertex->SetIndex(1,esdindex);
815 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
817 bestarray = new TObjArray(5);
818 fBestHypothesys.AddAt(bestarray,esdindex);
822 //setup tree of the prolongations
824 static AliITStrackMI tracks[7][100];
825 AliITStrackMI *currenttrack;
826 static AliITStrackMI currenttrack1;
827 static AliITStrackMI currenttrack2;
828 static AliITStrackMI backuptrack;
830 Int_t nindexes[7][100];
831 Float_t normalizedchi2[100];
832 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
833 otrack->SetNSkipped(0);
834 new (&(tracks[6][0])) AliITStrackMI(*otrack);
836 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
837 Int_t modstatus = 1; // found
841 // follow prolongations
842 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
845 AliITSlayer &layer=fgLayers[ilayer];
846 Double_t r = layer.GetR();
852 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
854 if (ntracks[ilayer]>=100) break;
855 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
856 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
857 if (ntracks[ilayer]>15+ilayer){
858 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
859 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
862 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
864 // material between SSD and SDD, SDD and SPD
866 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
868 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
872 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
873 Int_t idet=layer.FindDetectorIndex(phi,z);
875 Double_t trackGlobXYZ1[3];
876 currenttrack1.GetXYZ(trackGlobXYZ1);
878 // Get the budget to the primary vertex for the current track being prolonged
879 Double_t budgetToPrimVertex = GetEffectiveThickness();
881 // check if we allow a prolongation without point
882 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
884 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
885 // propagate to the layer radius
886 Double_t xToGo; vtrack->GetLocalXat(r,xToGo);
887 vtrack->AliExternalTrackParam::PropagateTo(xToGo,GetBz());
888 // apply correction for material of the current layer
889 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
890 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
891 vtrack->SetClIndex(ilayer,0);
892 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
893 LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc); // local module coords
894 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
895 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
900 // track outside layer acceptance in z
901 if (idet<0) continue;
903 //propagate to the intersection with the detector plane
904 const AliITSdetector &det=layer.GetDetector(idet);
905 new(¤ttrack2) AliITStrackMI(currenttrack1);
906 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
907 LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc); // local module coords
908 currenttrack2.Propagate(det.GetPhi(),det.GetR());
909 currenttrack1.SetDetectorIndex(idet);
910 currenttrack2.SetDetectorIndex(idet);
913 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
915 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
916 TMath::Sqrt(currenttrack1.GetSigmaZ2() +
917 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
918 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
919 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
920 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
921 TMath::Sqrt(currenttrack1.GetSigmaY2() +
922 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
923 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
924 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
926 // track at boundary between detectors, enlarge road
927 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
928 if ( (currenttrack1.GetY()-dy < det.GetYmin()+boundaryWidth) ||
929 (currenttrack1.GetY()+dy > det.GetYmax()-boundaryWidth) ||
930 (currenttrack1.GetZ()-dz < det.GetZmin()+boundaryWidth) ||
931 (currenttrack1.GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
932 Float_t tgl = TMath::Abs(currenttrack1.GetTgl());
933 if (tgl > 1.) tgl=1.;
934 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
935 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
936 Float_t snp = TMath::Abs(currenttrack1.GetSnp());
937 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) continue;
938 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
941 // road in global (rphi,z) [i.e. in tracking ref. system]
942 Double_t zmin = currenttrack1.GetZ() - dz;
943 Double_t zmax = currenttrack1.GetZ() + dz;
944 Double_t ymin = currenttrack1.GetY() + r*det.GetPhi() - dy;
945 Double_t ymax = currenttrack1.GetY() + r*det.GetPhi() + dy;
947 // select clusters in road
948 layer.SelectClusters(zmin,zmax,ymin,ymax);
949 //********************
951 // Define criteria for track-cluster association
952 Double_t msz = currenttrack1.GetSigmaZ2() +
953 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
954 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
955 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
956 Double_t msy = currenttrack1.GetSigmaY2() +
957 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
958 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
959 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
961 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
962 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
964 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
965 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
967 msz = 1./msz; // 1/RoadZ^2
968 msy = 1./msy; // 1/RoadY^2
971 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
973 const AliITSRecPoint *cl=0;
975 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
976 Bool_t deadzoneSPD=kFALSE;
977 currenttrack = ¤ttrack1;
979 // check if the road contains a dead zone
980 Int_t dead = CheckDeadZone(ilayer,idet,zmin,zmax);
981 // create a prolongation without clusters (check also if there are no clusters in the road)
983 ((layer.GetNextCluster(clidx,kTRUE))==0 &&
984 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
985 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
986 updatetrack->SetClIndex(ilayer,0);
988 modstatus = 5; // no cls in road
989 } else if (dead==1) {
990 modstatus = 7; // holes in z in SPD
991 } else if (dead==2) {
992 modstatus = 2; // dead from OCDB
994 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
995 // apply correction for material of the current layer
996 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
997 if (constrain) { // apply vertex constrain
998 updatetrack->SetConstrain(constrain);
999 Bool_t isPrim = kTRUE;
1000 if (ilayer<4) { // check that it's close to the vertex
1001 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1002 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1003 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1004 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1005 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1007 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1010 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1011 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1012 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1020 // loop over clusters in the road
1021 while ((cl=layer.GetNextCluster(clidx))!=0) {
1022 if (ntracks[ilayer]>95) break; //space for skipped clusters
1023 Bool_t changedet =kFALSE;
1024 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1025 Int_t idet=cl->GetDetectorIndex();
1027 if (currenttrack->GetDetectorIndex()==idet) { // track already on the cluster's detector
1028 // a first cut on track-cluster distance
1029 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1030 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1031 continue; // cluster not associated to track
1032 } else { // have to move track to cluster's detector
1033 const AliITSdetector &det=layer.GetDetector(idet);
1034 // a first cut on track-cluster distance
1036 if (!currenttrack2.GetProlongationFast(det.GetPhi(),det.GetR(),y,z)) continue;
1037 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1038 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1039 continue; // cluster not associated to track
1041 new (&backuptrack) AliITStrackMI(currenttrack2);
1043 currenttrack =¤ttrack2;
1044 if (!currenttrack->Propagate(det.GetPhi(),det.GetR())) {
1045 new (currenttrack) AliITStrackMI(backuptrack);
1049 currenttrack->SetDetectorIndex(idet);
1050 // Get again the budget to the primary vertex
1051 // for the current track being prolonged, if had to change detector
1052 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1055 // calculate track-clusters chi2
1056 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1058 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1059 if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1060 if (ntracks[ilayer]>=100) continue;
1061 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1062 updatetrack->SetClIndex(ilayer,0);
1063 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1065 if (cl->GetQ()!=0) { // real cluster
1066 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) continue;
1067 updatetrack->SetSampledEdx(cl->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
1068 modstatus = 1; // found
1069 } else { // virtual cluster in dead zone
1070 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1071 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1072 modstatus = 7; // holes in z in SPD
1076 Float_t xlocnewdet,zlocnewdet;
1077 LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet); // local module coords
1078 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1080 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1082 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1084 // apply correction for material of the current layer
1085 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1087 if (constrain) { // apply vertex constrain
1088 updatetrack->SetConstrain(constrain);
1089 Bool_t isPrim = kTRUE;
1090 if (ilayer<4) { // check that it's close to the vertex
1091 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1092 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1093 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1094 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1095 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1097 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1098 } //apply vertex constrain
1100 } // create new hypothesis
1101 } // loop over possible prolongations
1103 // allow one prolongation without clusters
1104 if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1105 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1106 // apply correction for material of the current layer
1107 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1108 vtrack->SetClIndex(ilayer,0);
1109 modstatus = 3; // skipped
1110 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1111 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1112 vtrack->IncrementNSkipped();
1116 // allow one prolongation without clusters for tracks with |tgl|>1.1
1117 if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) { //big theta - for low flux
1118 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1119 // apply correction for material of the current layer
1120 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1121 vtrack->SetClIndex(ilayer,0);
1122 modstatus = 3; // skipped
1123 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1124 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1125 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1130 } // loop over tracks in layer ilayer+1
1132 //loop over track candidates for the current layer
1138 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1139 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1140 if (normalizedchi2[itrack] <
1141 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1145 if (constrain) { // constrain
1146 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1148 } else { // !constrain
1149 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1154 // sort tracks by increasing normalized chi2
1155 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1156 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1157 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1158 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1159 } // end loop over layers
1163 // Now select tracks to be kept
1165 Int_t max = constrain ? 20 : 5;
1167 // tracks that reach layer 0 (SPD inner)
1168 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1169 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1170 if (track.GetNumberOfClusters()<2) continue;
1171 if (!constrain && track.GetNormChi2(0) >
1172 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1175 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1178 // tracks that reach layer 1 (SPD outer)
1179 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1180 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1181 if (track.GetNumberOfClusters()<4) continue;
1182 if (!constrain && track.GetNormChi2(1) >
1183 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1184 if (constrain) track.IncrementNSkipped();
1186 track.SetD(0,track.GetD(GetX(),GetY()));
1187 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1188 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1189 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1192 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1195 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1197 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1198 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1199 if (track.GetNumberOfClusters()<3) continue;
1200 if (!constrain && track.GetNormChi2(2) >
1201 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1202 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1204 track.SetD(0,track.GetD(GetX(),GetY()));
1205 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1206 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1207 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1210 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1216 // register best track of each layer - important for V0 finder
1218 for (Int_t ilayer=0;ilayer<5;ilayer++){
1219 if (ntracks[ilayer]==0) continue;
1220 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1221 if (track.GetNumberOfClusters()<1) continue;
1222 CookLabel(&track,0);
1223 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1227 // update TPC V0 information
1229 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1230 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1231 for (Int_t i=0;i<3;i++){
1232 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1233 if (index==0) break;
1234 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1235 if (vertex->GetStatus()<0) continue; // rejected V0
1237 if (otrack->GetSign()>0) {
1238 vertex->SetIndex(0,esdindex);
1241 vertex->SetIndex(1,esdindex);
1243 //find nearest layer with track info
1244 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1245 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1246 Int_t nearest = nearestold;
1247 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1248 if (ntracks[nearest]==0){
1253 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1254 if (nearestold<5&&nearest<5){
1255 Bool_t accept = track.GetNormChi2(nearest)<10;
1257 if (track.GetSign()>0) {
1258 vertex->SetParamP(track);
1259 vertex->Update(fprimvertex);
1260 //vertex->SetIndex(0,track.fESDtrack->GetID());
1261 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1263 vertex->SetParamN(track);
1264 vertex->Update(fprimvertex);
1265 //vertex->SetIndex(1,track.fESDtrack->GetID());
1266 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1268 vertex->SetStatus(vertex->GetStatus()+1);
1270 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1277 //------------------------------------------------------------------------
1278 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1280 //--------------------------------------------------------------------
1283 return fgLayers[layer];
1285 //------------------------------------------------------------------------
1286 AliITStrackerMI::AliITSlayer::AliITSlayer():
1311 //--------------------------------------------------------------------
1312 //default AliITSlayer constructor
1313 //--------------------------------------------------------------------
1314 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1315 fClusterWeight[i]=0;
1316 fClusterTracks[0][i]=-1;
1317 fClusterTracks[1][i]=-1;
1318 fClusterTracks[2][i]=-1;
1319 fClusterTracks[3][i]=-1;
1322 //------------------------------------------------------------------------
1323 AliITStrackerMI::AliITSlayer::
1324 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1349 //--------------------------------------------------------------------
1350 //main AliITSlayer constructor
1351 //--------------------------------------------------------------------
1352 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1353 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1355 //------------------------------------------------------------------------
1356 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1358 fPhiOffset(layer.fPhiOffset),
1359 fNladders(layer.fNladders),
1360 fZOffset(layer.fZOffset),
1361 fNdetectors(layer.fNdetectors),
1362 fDetectors(layer.fDetectors),
1367 fClustersCs(layer.fClustersCs),
1368 fClusterIndexCs(layer.fClusterIndexCs),
1372 fCurrentSlice(layer.fCurrentSlice),
1379 fAccepted(layer.fAccepted),
1383 //------------------------------------------------------------------------
1384 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1385 //--------------------------------------------------------------------
1386 // AliITSlayer destructor
1387 //--------------------------------------------------------------------
1388 delete[] fDetectors;
1389 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1390 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1391 fClusterWeight[i]=0;
1392 fClusterTracks[0][i]=-1;
1393 fClusterTracks[1][i]=-1;
1394 fClusterTracks[2][i]=-1;
1395 fClusterTracks[3][i]=-1;
1398 //------------------------------------------------------------------------
1399 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1400 //--------------------------------------------------------------------
1401 // This function removes loaded clusters
1402 //--------------------------------------------------------------------
1403 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1404 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1405 fClusterWeight[i]=0;
1406 fClusterTracks[0][i]=-1;
1407 fClusterTracks[1][i]=-1;
1408 fClusterTracks[2][i]=-1;
1409 fClusterTracks[3][i]=-1;
1415 //------------------------------------------------------------------------
1416 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1417 //--------------------------------------------------------------------
1418 // This function reset weights of the clusters
1419 //--------------------------------------------------------------------
1420 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1421 fClusterWeight[i]=0;
1422 fClusterTracks[0][i]=-1;
1423 fClusterTracks[1][i]=-1;
1424 fClusterTracks[2][i]=-1;
1425 fClusterTracks[3][i]=-1;
1427 for (Int_t i=0; i<fN;i++) {
1428 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1429 if (cl&&cl->IsUsed()) cl->Use();
1433 //------------------------------------------------------------------------
1434 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1435 //--------------------------------------------------------------------
1436 // This function calculates the road defined by the cluster density
1437 //--------------------------------------------------------------------
1439 for (Int_t i=0; i<fN; i++) {
1440 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1442 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1444 //------------------------------------------------------------------------
1445 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1446 //--------------------------------------------------------------------
1447 //This function adds a cluster to this layer
1448 //--------------------------------------------------------------------
1449 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1450 ::Error("InsertCluster","Too many clusters !\n");
1456 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1457 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1458 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1459 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1460 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1464 //------------------------------------------------------------------------
1465 void AliITStrackerMI::AliITSlayer::SortClusters()
1470 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1471 Float_t *z = new Float_t[fN];
1472 Int_t * index = new Int_t[fN];
1474 for (Int_t i=0;i<fN;i++){
1475 z[i] = fClusters[i]->GetZ();
1477 TMath::Sort(fN,z,index,kFALSE);
1478 for (Int_t i=0;i<fN;i++){
1479 clusters[i] = fClusters[index[i]];
1482 for (Int_t i=0;i<fN;i++){
1483 fClusters[i] = clusters[i];
1484 fZ[i] = fClusters[i]->GetZ();
1485 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1486 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1487 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1497 for (Int_t i=0;i<fN;i++){
1498 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1499 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1500 fClusterIndex[i] = i;
1504 fDy5 = (fYB[1]-fYB[0])/5.;
1505 fDy10 = (fYB[1]-fYB[0])/10.;
1506 fDy20 = (fYB[1]-fYB[0])/20.;
1507 for (Int_t i=0;i<6;i++) fN5[i] =0;
1508 for (Int_t i=0;i<11;i++) fN10[i]=0;
1509 for (Int_t i=0;i<21;i++) fN20[i]=0;
1511 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;}
1512 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;}
1513 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;}
1516 for (Int_t i=0;i<fN;i++)
1517 for (Int_t irot=-1;irot<=1;irot++){
1518 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1520 for (Int_t slice=0; slice<6;slice++){
1521 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1522 fClusters5[slice][fN5[slice]] = fClusters[i];
1523 fY5[slice][fN5[slice]] = curY;
1524 fZ5[slice][fN5[slice]] = fZ[i];
1525 fClusterIndex5[slice][fN5[slice]]=i;
1530 for (Int_t slice=0; slice<11;slice++){
1531 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1532 fClusters10[slice][fN10[slice]] = fClusters[i];
1533 fY10[slice][fN10[slice]] = curY;
1534 fZ10[slice][fN10[slice]] = fZ[i];
1535 fClusterIndex10[slice][fN10[slice]]=i;
1540 for (Int_t slice=0; slice<21;slice++){
1541 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1542 fClusters20[slice][fN20[slice]] = fClusters[i];
1543 fY20[slice][fN20[slice]] = curY;
1544 fZ20[slice][fN20[slice]] = fZ[i];
1545 fClusterIndex20[slice][fN20[slice]]=i;
1552 // consistency check
1554 for (Int_t i=0;i<fN-1;i++){
1560 for (Int_t slice=0;slice<21;slice++)
1561 for (Int_t i=0;i<fN20[slice]-1;i++){
1562 if (fZ20[slice][i]>fZ20[slice][i+1]){
1569 //------------------------------------------------------------------------
1570 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1571 //--------------------------------------------------------------------
1572 // This function returns the index of the nearest cluster
1573 //--------------------------------------------------------------------
1576 if (fCurrentSlice<0) {
1585 if (ncl==0) return 0;
1586 Int_t b=0, e=ncl-1, m=(b+e)/2;
1587 for (; b<e; m=(b+e)/2) {
1588 // if (z > fClusters[m]->GetZ()) b=m+1;
1589 if (z > zcl[m]) b=m+1;
1594 //------------------------------------------------------------------------
1595 void AliITStrackerMI::AliITSlayer::
1596 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1597 //--------------------------------------------------------------------
1598 // This function sets the "window"
1599 //--------------------------------------------------------------------
1601 Double_t circle=2*TMath::Pi()*fR;
1602 fYmin = ymin; fYmax =ymax;
1603 Float_t ymiddle = (fYmin+fYmax)*0.5;
1604 if (ymiddle<fYB[0]) {
1605 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1606 } else if (ymiddle>fYB[1]) {
1607 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1613 fClustersCs = fClusters;
1614 fClusterIndexCs = fClusterIndex;
1620 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1621 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1622 if (slice<0) slice=0;
1623 if (slice>20) slice=20;
1624 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1626 fCurrentSlice=slice;
1627 fClustersCs = fClusters20[fCurrentSlice];
1628 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1629 fYcs = fY20[fCurrentSlice];
1630 fZcs = fZ20[fCurrentSlice];
1631 fNcs = fN20[fCurrentSlice];
1636 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1637 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1638 if (slice<0) slice=0;
1639 if (slice>10) slice=10;
1640 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1642 fCurrentSlice=slice;
1643 fClustersCs = fClusters10[fCurrentSlice];
1644 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1645 fYcs = fY10[fCurrentSlice];
1646 fZcs = fZ10[fCurrentSlice];
1647 fNcs = fN10[fCurrentSlice];
1652 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1653 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1654 if (slice<0) slice=0;
1655 if (slice>5) slice=5;
1656 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1658 fCurrentSlice=slice;
1659 fClustersCs = fClusters5[fCurrentSlice];
1660 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1661 fYcs = fY5[fCurrentSlice];
1662 fZcs = fZ5[fCurrentSlice];
1663 fNcs = fN5[fCurrentSlice];
1667 fI=FindClusterIndex(zmin); fZmax=zmax;
1668 fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1674 //------------------------------------------------------------------------
1675 Int_t AliITStrackerMI::AliITSlayer::
1676 FindDetectorIndex(Double_t phi, Double_t z) const {
1677 //--------------------------------------------------------------------
1678 //This function finds the detector crossed by the track
1679 //--------------------------------------------------------------------
1681 if (fZOffset<0) // old geometry
1682 dphi = -(phi-fPhiOffset);
1683 else // new geometry
1684 dphi = phi-fPhiOffset;
1687 if (dphi < 0) dphi += 2*TMath::Pi();
1688 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1689 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1690 if (np>=fNladders) np-=fNladders;
1691 if (np<0) np+=fNladders;
1694 Double_t dz=fZOffset-z;
1695 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1696 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1697 if (nz>=fNdetectors) return -1;
1698 if (nz<0) return -1;
1700 return np*fNdetectors + nz;
1702 //------------------------------------------------------------------------
1703 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1705 //--------------------------------------------------------------------
1706 // This function returns clusters within the "window"
1707 //--------------------------------------------------------------------
1709 if (fCurrentSlice<0) {
1710 Double_t rpi2 = 2.*fR*TMath::Pi();
1711 for (Int_t i=fI; i<fImax; i++) {
1713 if (fYmax<y) y -= rpi2;
1714 if (fYmin>y) y += rpi2;
1715 if (y<fYmin) continue;
1716 if (y>fYmax) continue;
1717 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1720 return fClusters[i];
1723 for (Int_t i=fI; i<fImax; i++) {
1724 if (fYcs[i]<fYmin) continue;
1725 if (fYcs[i]>fYmax) continue;
1726 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1727 ci=fClusterIndexCs[i];
1729 return fClustersCs[i];
1734 //------------------------------------------------------------------------
1735 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1737 //--------------------------------------------------------------------
1738 // This function returns the layer thickness at this point (units X0)
1739 //--------------------------------------------------------------------
1741 x0=AliITSRecoParam::GetX0Air();
1742 if (43<fR&&fR<45) { //SSD2
1745 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1746 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1747 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1748 for (Int_t i=0; i<12; i++) {
1749 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1750 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1754 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1755 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1759 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1760 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1763 if (37<fR&&fR<41) { //SSD1
1766 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1767 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1768 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1769 for (Int_t i=0; i<11; i++) {
1770 if (TMath::Abs(z-3.9*i)<0.15) {
1771 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1775 if (TMath::Abs(z+3.9*i)<0.15) {
1776 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1780 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1781 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1784 if (13<fR&&fR<26) { //SDD
1787 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1789 if (TMath::Abs(y-1.80)<0.55) {
1791 for (Int_t j=0; j<20; j++) {
1792 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1793 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1796 if (TMath::Abs(y+1.80)<0.55) {
1798 for (Int_t j=0; j<20; j++) {
1799 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1800 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1804 for (Int_t i=0; i<4; i++) {
1805 if (TMath::Abs(z-7.3*i)<0.60) {
1807 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1810 if (TMath::Abs(z+7.3*i)<0.60) {
1812 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1817 if (6<fR&&fR<8) { //SPD2
1818 Double_t dd=0.0063; x0=21.5;
1820 if (TMath::Abs(y-3.08)>0.5) d+=dd;
1821 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
1823 if (3<fR&&fR<5) { //SPD1
1824 Double_t dd=0.0063; x0=21.5;
1826 if (TMath::Abs(y+0.21)>0.6) d+=dd;
1827 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
1832 //------------------------------------------------------------------------
1833 Double_t AliITStrackerMI::GetEffectiveThickness()
1835 //--------------------------------------------------------------------
1836 // Returns the thickness between the current layer and the vertex (units X0)
1837 //--------------------------------------------------------------------
1840 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
1841 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
1842 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
1846 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
1847 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
1851 Double_t xn=fgLayers[fI].GetR();
1852 for (Int_t i=0; i<fI; i++) {
1853 Double_t xi=fgLayers[i].GetR();
1854 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
1860 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
1861 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
1864 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
1865 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
1869 //------------------------------------------------------------------------
1870 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
1871 //-------------------------------------------------------------------
1872 // This function returns number of clusters within the "window"
1873 //--------------------------------------------------------------------
1875 for (Int_t i=fI; i<fN; i++) {
1876 const AliITSRecPoint *c=fClusters[i];
1877 if (c->GetZ() > fZmax) break;
1878 if (c->IsUsed()) continue;
1879 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1880 Double_t y=fR*det.GetPhi() + c->GetY();
1882 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
1883 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1885 if (y<fYmin) continue;
1886 if (y>fYmax) continue;
1891 //------------------------------------------------------------------------
1892 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
1893 const AliITStrackMI *clusters,Bool_t extra)
1895 //--------------------------------------------------------------------
1896 // This function refits the track "track" at the position "x" using
1897 // the clusters from "clusters"
1898 // If "extra"==kTRUE,
1899 // the clusters from overlapped modules get attached to "track"
1900 //--------------------------------------------------------------------
1902 Int_t index[AliITSgeomTGeo::kNLayers];
1904 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
1905 Int_t nc=clusters->GetNumberOfClusters();
1906 for (k=0; k<nc; k++) {
1907 Int_t idx=clusters->GetClusterIndex(k);
1908 Int_t ilayer=(idx&0xf0000000)>>28;
1912 return RefitAt(xx,track,index,extra); // call the method below
1914 //------------------------------------------------------------------------
1915 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
1916 const Int_t *clusters,Bool_t extra)
1918 //--------------------------------------------------------------------
1919 // This function refits the track "track" at the position "x" using
1920 // the clusters from array
1921 // If "extra"==kTRUE,
1922 // the clusters from overlapped modules get attached to "track"
1923 //--------------------------------------------------------------------
1924 Int_t index[AliITSgeomTGeo::kNLayers];
1926 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
1928 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
1929 index[k]=clusters[k];
1932 // special for cosmics: check which the innermost layer crossed
1934 Int_t innermostlayer=5;
1935 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
1936 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
1937 if(drphi < fgLayers[innermostlayer].GetR()) break;
1939 //printf(" drphi %f innermost %d\n",drphi,innermostlayer);
1941 Int_t modstatus=1; // found
1943 Int_t from, to, step;
1944 if (xx > track->GetX()) {
1945 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
1948 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
1951 TString dir = (step>0 ? "outward" : "inward");
1953 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
1954 AliITSlayer &layer=fgLayers[ilayer];
1955 Double_t r=layer.GetR();
1956 if (step<0 && xx>r) break;
1958 // material between SSD and SDD, SDD and SPD
1959 Double_t hI=ilayer-0.5*step;
1960 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
1961 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
1962 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
1963 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
1965 // remember old position [SR, GSI 18.02.2003]
1966 Double_t oldX=0., oldY=0., oldZ=0.;
1967 if (track->IsStartedTimeIntegral() && step==1) {
1968 track->GetGlobalXYZat(track->GetX(),oldX,oldY,oldZ);
1972 Double_t oldGlobXYZ[3];
1973 track->GetXYZ(oldGlobXYZ);
1976 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
1978 Int_t idet=layer.FindDetectorIndex(phi,z);
1980 // check if we allow a prolongation without point for large-eta tracks
1981 Int_t skip = CheckSkipLayer(track,ilayer,idet);
1983 // propagate to the layer radius
1984 Double_t xToGo; track->GetLocalXat(r,xToGo);
1985 track->AliExternalTrackParam::PropagateTo(xToGo,GetBz());
1986 // apply correction for material of the current layer
1987 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
1988 modstatus = 4; // out in z
1989 LocalModuleCoord(ilayer,idet,track,xloc,zloc); // local module coords
1990 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1991 // track time update [SR, GSI 17.02.2003]
1992 if (track->IsStartedTimeIntegral() && step==1) {
1993 Double_t newX, newY, newZ;
1994 track->GetGlobalXYZat(track->GetX(),newX,newY,newZ);
1995 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1996 (oldZ-newZ)*(oldZ-newZ);
1997 track->AddTimeStep(TMath::Sqrt(dL2));
2002 if (idet<0) return kFALSE;
2004 const AliITSdetector &det=layer.GetDetector(idet);
2006 if (!track->Propagate(phi,det.GetR())) return kFALSE;
2007 track->SetDetectorIndex(idet);
2008 LocalModuleCoord(ilayer,idet,track,xloc,zloc); // local module coords
2010 Double_t dz,zmin,zmax;
2012 const AliITSRecPoint *clAcc=0;
2013 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2015 Int_t idx=index[ilayer];
2016 if (idx>=0) { // cluster in this layer
2017 modstatus = 6; // no refit
2018 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2020 if (idet != cl->GetDetectorIndex()) {
2021 idet=cl->GetDetectorIndex();
2022 const AliITSdetector &det=layer.GetDetector(idet);
2023 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2024 track->SetDetectorIndex(idet);
2025 LocalModuleCoord(ilayer,idet,track,xloc,zloc); // local module coords
2027 //Double_t chi2=track->GetPredictedChi2(cl);
2028 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2029 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2033 modstatus = 1; // found
2038 } else { // no cluster in this layer
2040 modstatus = 3; // skipped
2042 modstatus = 5; // no cls in road
2044 dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
2045 TMath::Sqrt(track->GetSigmaZ2() +
2046 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
2047 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
2048 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
2049 zmin=track->GetZ() - dz;
2050 zmax=track->GetZ() + dz;
2051 Int_t dead = CheckDeadZone(ilayer,idet,zmin,zmax);
2052 if (dead==1) modstatus = 7; // holes in z in SPD
2053 if (dead==2) modstatus = 2; // dead from OCDB
2058 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2059 track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2061 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2064 if (extra) { // search for extra clusters in overlapped modules
2065 AliITStrackV2 tmp(*track);
2066 Double_t dy,ymin,ymax;
2067 dz=4*TMath::Sqrt(tmp.GetSigmaZ2()+AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
2068 if (dz < 0.5*TMath::Abs(tmp.GetTgl())) dz=0.5*TMath::Abs(tmp.GetTgl());
2069 dy=4*TMath::Sqrt(track->GetSigmaY2()+AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
2070 if (dy < 0.5*TMath::Abs(tmp.GetSnp())) dy=0.5*TMath::Abs(tmp.GetSnp());
2071 zmin=track->GetZ() - dz;
2072 zmax=track->GetZ() + dz;
2073 ymin=track->GetY() + phi*r - dy;
2074 ymax=track->GetY() + phi*r + dy;
2075 layer.SelectClusters(zmin,zmax,ymin,ymax);
2077 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2079 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2(), tolerance=0.1;
2080 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2081 // only clusters in another module! (overlaps)
2082 idetExtra = clExtra->GetDetectorIndex();
2083 if (idet == idetExtra) continue;
2085 const AliITSdetector &det=layer.GetDetector(idetExtra);
2087 if (!tmp.Propagate(det.GetPhi(),det.GetR())) continue;
2089 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2090 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2092 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2093 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2096 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2097 track->SetExtraModule(ilayer,idetExtra);
2099 } // end search for extra clusters in overlapped modules
2101 // Correct for material of the current layer
2102 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2104 // track time update [SR, GSI 17.02.2003]
2105 if (track->IsStartedTimeIntegral() && step==1) {
2106 Double_t newX, newY, newZ;
2107 track->GetGlobalXYZat(track->GetX(),newX,newY,newZ);
2108 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2109 (oldZ-newZ)*(oldZ-newZ);
2110 track->AddTimeStep(TMath::Sqrt(dL2));
2114 } // end loop on layers
2116 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2120 //------------------------------------------------------------------------
2121 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2124 // calculate normalized chi2
2125 // return NormalizedChi2(track,0);
2128 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2129 // track->fdEdxMismatch=0;
2130 Float_t dedxmismatch =0;
2131 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2133 for (Int_t i = 0;i<6;i++){
2134 if (track->GetClIndex(i)>0){
2135 Float_t cerry, cerrz;
2136 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2138 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2141 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2142 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2143 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2145 cchi2+=(0.5-ratio)*10.;
2146 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2147 dedxmismatch+=(0.5-ratio)*10.;
2151 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2152 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2153 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2154 if (i<2) chi2+=2*cl->GetDeltaProbability();
2160 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2161 track->SetdEdxMismatch(dedxmismatch);
2165 for (Int_t i = 0;i<4;i++){
2166 if (track->GetClIndex(i)>0){
2167 Float_t cerry, cerrz;
2168 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2169 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2172 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2173 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2177 for (Int_t i = 4;i<6;i++){
2178 if (track->GetClIndex(i)>0){
2179 Float_t cerry, cerrz;
2180 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2181 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2184 Float_t cerryb, cerrzb;
2185 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2186 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2189 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2190 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2195 if (track->GetESDtrack()->GetTPCsignal()>85){
2196 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2198 chi2+=(0.5-ratio)*5.;
2201 chi2+=(ratio-2.0)*3;
2205 Double_t match = TMath::Sqrt(track->GetChi22());
2206 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2207 if (!track->GetConstrain()) {
2208 if (track->GetNumberOfClusters()>2) {
2209 match/=track->GetNumberOfClusters()-2.;
2214 if (match<0) match=0;
2215 Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2216 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2217 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2218 1./(1.+track->GetNSkipped()));
2222 //------------------------------------------------------------------------
2223 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2226 // return matching chi2 between two tracks
2227 AliITStrackMI track3(*track2);
2228 track3.Propagate(track1->GetAlpha(),track1->GetX());
2230 vec(0,0)=track1->GetY() - track3.GetY();
2231 vec(1,0)=track1->GetZ() - track3.GetZ();
2232 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2233 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2234 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2237 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2238 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2239 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2240 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2241 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2243 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2244 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2245 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2246 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2248 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2249 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2250 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2252 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2253 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2255 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2258 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2259 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2262 //------------------------------------------------------------------------
2263 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr)
2266 // return probability that given point (characterized by z position and error)
2267 // is in SPD dead zone
2269 Double_t probability = 0.;
2270 Double_t absz = TMath::Abs(zpos);
2271 Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2272 0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2273 if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2274 Double_t zmin, zmax;
2275 if (zpos<-6.) { // dead zone at z = -7
2276 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2277 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2278 } else if (zpos>6.) { // dead zone at z = +7
2279 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2280 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2281 } else if (absz<2.) { // dead zone at z = 0
2282 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2283 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2288 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2290 probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2291 TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2294 //------------------------------------------------------------------------
2295 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2298 // calculate normalized chi2
2300 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2302 for (Int_t i = 0;i<6;i++){
2303 if (TMath::Abs(track->GetDy(i))>0){
2304 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2305 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2308 else{chi2[i]=10000;}
2311 TMath::Sort(6,chi2,index,kFALSE);
2312 Float_t max = float(ncl)*fac-1.;
2313 Float_t sumchi=0, sumweight=0;
2314 for (Int_t i=0;i<max+1;i++){
2315 Float_t weight = (i<max)?1.:(max+1.-i);
2316 sumchi+=weight*chi2[index[i]];
2319 Double_t normchi2 = sumchi/sumweight;
2322 //------------------------------------------------------------------------
2323 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2326 // calculate normalized chi2
2327 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2330 for (Int_t i=0;i<6;i++){
2331 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2332 Double_t sy1 = forwardtrack->GetSigmaY(i);
2333 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2334 Double_t sy2 = backtrack->GetSigmaY(i);
2335 Double_t sz2 = backtrack->GetSigmaZ(i);
2336 if (i<2){ sy2=1000.;sz2=1000;}
2338 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2339 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2341 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2342 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2344 res+= nz0*nz0+ny0*ny0;
2347 if (npoints>1) return
2348 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2349 //2*forwardtrack->fNUsed+
2350 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2351 1./(1.+forwardtrack->GetNSkipped()));
2354 //------------------------------------------------------------------------
2355 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2356 //--------------------------------------------------------------------
2357 // Return pointer to a given cluster
2358 //--------------------------------------------------------------------
2359 Int_t l=(index & 0xf0000000) >> 28;
2360 Int_t c=(index & 0x0fffffff) >> 00;
2361 return fgLayers[l].GetWeight(c);
2363 //------------------------------------------------------------------------
2364 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2366 //---------------------------------------------
2367 // register track to the list
2369 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2372 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2373 Int_t index = track->GetClusterIndex(icluster);
2374 Int_t l=(index & 0xf0000000) >> 28;
2375 Int_t c=(index & 0x0fffffff) >> 00;
2376 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2377 for (Int_t itrack=0;itrack<4;itrack++){
2378 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2379 fgLayers[l].SetClusterTracks(itrack,c,id);
2385 //------------------------------------------------------------------------
2386 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2388 //---------------------------------------------
2389 // unregister track from the list
2390 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2391 Int_t index = track->GetClusterIndex(icluster);
2392 Int_t l=(index & 0xf0000000) >> 28;
2393 Int_t c=(index & 0x0fffffff) >> 00;
2394 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2395 for (Int_t itrack=0;itrack<4;itrack++){
2396 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2397 fgLayers[l].SetClusterTracks(itrack,c,-1);
2402 //------------------------------------------------------------------------
2403 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2405 //-------------------------------------------------------------
2406 //get number of shared clusters
2407 //-------------------------------------------------------------
2409 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2410 // mean number of clusters
2411 Float_t *ny = GetNy(id), *nz = GetNz(id);
2414 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2415 Int_t index = track->GetClusterIndex(icluster);
2416 Int_t l=(index & 0xf0000000) >> 28;
2417 Int_t c=(index & 0x0fffffff) >> 00;
2418 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2420 printf("problem\n");
2422 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2426 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2427 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2428 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2430 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2433 deltan = (cl->GetNz()-nz[l]);
2435 if (deltan>2.0) continue; // extended - highly probable shared cluster
2436 weight = 2./TMath::Max(3.+deltan,2.);
2438 for (Int_t itrack=0;itrack<4;itrack++){
2439 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2441 clist[l] = (AliITSRecPoint*)GetCluster(index);
2447 track->SetNUsed(shared);
2450 //------------------------------------------------------------------------
2451 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2454 // find first shared track
2456 // mean number of clusters
2457 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2459 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2460 Int_t sharedtrack=100000;
2461 Int_t tracks[24],trackindex=0;
2462 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2464 for (Int_t icluster=0;icluster<6;icluster++){
2465 if (clusterlist[icluster]<0) continue;
2466 Int_t index = clusterlist[icluster];
2467 Int_t l=(index & 0xf0000000) >> 28;
2468 Int_t c=(index & 0x0fffffff) >> 00;
2470 printf("problem\n");
2472 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2473 //if (l>3) continue;
2474 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2477 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2478 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2479 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2481 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2484 deltan = (cl->GetNz()-nz[l]);
2486 if (deltan>2.0) continue; // extended - highly probable shared cluster
2488 for (Int_t itrack=3;itrack>=0;itrack--){
2489 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2490 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2491 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2496 if (trackindex==0) return -1;
2498 sharedtrack = tracks[0];
2500 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2503 Int_t track[24], cluster[24];
2504 for (Int_t i=0;i<trackindex;i++){ track[i]=-1; cluster[i]=0;}
2507 for (Int_t i=0;i<trackindex;i++){
2508 if (tracks[i]<0) continue;
2509 track[index] = tracks[i];
2511 for (Int_t j=i+1;j<trackindex;j++){
2512 if (tracks[j]<0) continue;
2513 if (tracks[j]==tracks[i]){
2521 for (Int_t i=0;i<index;i++){
2522 if (cluster[index]>max) {
2523 sharedtrack=track[index];
2530 if (sharedtrack>=100000) return -1;
2532 // make list of overlaps
2534 for (Int_t icluster=0;icluster<6;icluster++){
2535 if (clusterlist[icluster]<0) continue;
2536 Int_t index = clusterlist[icluster];
2537 Int_t l=(index & 0xf0000000) >> 28;
2538 Int_t c=(index & 0x0fffffff) >> 00;
2539 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2540 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2542 if (cl->GetNy()>2) continue;
2543 if (cl->GetNz()>2) continue;
2546 if (cl->GetNy()>3) continue;
2547 if (cl->GetNz()>3) continue;
2550 for (Int_t itrack=3;itrack>=0;itrack--){
2551 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2552 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2560 //------------------------------------------------------------------------
2561 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2563 // try to find track hypothesys without conflicts
2564 // with minimal chi2;
2565 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2566 Int_t entries1 = arr1->GetEntriesFast();
2567 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2568 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2569 Int_t entries2 = arr2->GetEntriesFast();
2570 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2572 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2573 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2574 if (track10->Pt()>0.5+track20->Pt()) return track10;
2576 for (Int_t itrack=0;itrack<entries1;itrack++){
2577 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2578 UnRegisterClusterTracks(track,trackID1);
2581 for (Int_t itrack=0;itrack<entries2;itrack++){
2582 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2583 UnRegisterClusterTracks(track,trackID2);
2587 Float_t maxconflicts=6;
2588 Double_t maxchi2 =1000.;
2590 // get weight of hypothesys - using best hypothesy
2593 Int_t list1[6],list2[6];
2594 AliITSRecPoint *clist1[6], *clist2[6] ;
2595 RegisterClusterTracks(track10,trackID1);
2596 RegisterClusterTracks(track20,trackID2);
2597 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2598 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2599 UnRegisterClusterTracks(track10,trackID1);
2600 UnRegisterClusterTracks(track20,trackID2);
2603 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2604 Float_t nerry[6],nerrz[6];
2605 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2606 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2607 for (Int_t i=0;i<6;i++){
2608 if ( (erry1[i]>0) && (erry2[i]>0)) {
2609 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2610 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2612 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2613 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2615 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2616 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2617 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2620 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2621 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2622 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2630 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2631 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2632 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2633 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2635 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2636 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2637 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2639 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2640 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2641 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2644 Double_t sumw = w1+w2;
2648 w1 = (d2+0.5)/(d1+d2+1.);
2649 w2 = (d1+0.5)/(d1+d2+1.);
2651 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2652 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2654 // get pair of "best" hypothesys
2656 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2657 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2659 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2660 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2661 //if (track1->fFakeRatio>0) continue;
2662 RegisterClusterTracks(track1,trackID1);
2663 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2664 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2666 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2667 //if (track2->fFakeRatio>0) continue;
2669 RegisterClusterTracks(track2,trackID2);
2670 Int_t list1[6],list2[6];
2671 AliITSRecPoint *clist1[6], *clist2[6] ;
2672 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2673 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2674 UnRegisterClusterTracks(track2,trackID2);
2676 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2677 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2678 if (nskipped>0.5) continue;
2680 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2681 if (conflict1+1<cconflict1) continue;
2682 if (conflict2+1<cconflict2) continue;
2686 for (Int_t i=0;i<6;i++){
2688 Float_t c1 =0.; // conflict coeficients
2690 if (clist1[i]&&clist2[i]){
2693 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2696 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2698 c1 = 2./TMath::Max(3.+deltan,2.);
2699 c2 = 2./TMath::Max(3.+deltan,2.);
2705 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2708 deltan = (clist1[i]->GetNz()-nz1[i]);
2710 c1 = 2./TMath::Max(3.+deltan,2.);
2717 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2720 deltan = (clist2[i]->GetNz()-nz2[i]);
2722 c2 = 2./TMath::Max(3.+deltan,2.);
2727 Double_t chi21=0,chi22=0;
2728 if (TMath::Abs(track1->GetDy(i))>0.) {
2729 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2730 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2731 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2732 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2734 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2737 if (TMath::Abs(track2->GetDy(i))>0.) {
2738 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2739 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2740 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2741 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2744 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2746 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2747 if (chi21>0) sum+=w1;
2748 if (chi22>0) sum+=w2;
2751 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2752 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2753 Double_t normchi2 = 2*conflict+sumchi2/sum;
2754 if ( normchi2 <maxchi2 ){
2757 maxconflicts = conflict;
2761 UnRegisterClusterTracks(track1,trackID1);
2764 // if (maxconflicts<4 && maxchi2<th0){
2765 if (maxchi2<th0*2.){
2766 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
2767 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2768 track1->SetChi2MIP(5,maxconflicts);
2769 track1->SetChi2MIP(6,maxchi2);
2770 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
2771 // track1->UpdateESDtrack(AliESDtrack::kITSin);
2772 track1->SetChi2MIP(8,index1);
2773 fBestTrackIndex[trackID1] =index1;
2774 UpdateESDtrack(track1, AliESDtrack::kITSin);
2776 else if (track10->GetChi2MIP(0)<th1){
2777 track10->SetChi2MIP(5,maxconflicts);
2778 track10->SetChi2MIP(6,maxchi2);
2779 // track10->UpdateESDtrack(AliESDtrack::kITSin);
2780 UpdateESDtrack(track10,AliESDtrack::kITSin);
2783 for (Int_t itrack=0;itrack<entries1;itrack++){
2784 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2785 UnRegisterClusterTracks(track,trackID1);
2788 for (Int_t itrack=0;itrack<entries2;itrack++){
2789 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2790 UnRegisterClusterTracks(track,trackID2);
2793 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
2794 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
2795 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
2796 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
2797 RegisterClusterTracks(track10,trackID1);
2799 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
2800 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
2801 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
2802 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
2803 RegisterClusterTracks(track20,trackID2);
2808 //------------------------------------------------------------------------
2809 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
2810 //--------------------------------------------------------------------
2811 // This function marks clusters assigned to the track
2812 //--------------------------------------------------------------------
2813 AliTracker::UseClusters(t,from);
2815 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
2816 //if (c->GetQ()>2) c->Use();
2817 if (c->GetSigmaZ2()>0.1) c->Use();
2818 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
2819 //if (c->GetQ()>2) c->Use();
2820 if (c->GetSigmaZ2()>0.1) c->Use();
2823 //------------------------------------------------------------------------
2824 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
2826 //------------------------------------------------------------------
2827 // add track to the list of hypothesys
2828 //------------------------------------------------------------------
2830 if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
2832 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2834 array = new TObjArray(10);
2835 fTrackHypothesys.AddAt(array,esdindex);
2837 array->AddLast(track);
2839 //------------------------------------------------------------------------
2840 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
2842 //-------------------------------------------------------------------
2843 // compress array of track hypothesys
2844 // keep only maxsize best hypothesys
2845 //-------------------------------------------------------------------
2846 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
2847 if (! (fTrackHypothesys.At(esdindex)) ) return;
2848 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2849 Int_t entries = array->GetEntriesFast();
2851 //- find preliminary besttrack as a reference
2852 Float_t minchi2=10000;
2854 AliITStrackMI * besttrack=0;
2855 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
2856 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2857 if (!track) continue;
2858 Float_t chi2 = NormalizedChi2(track,0);
2860 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
2861 track->SetLabel(tpcLabel);
2863 track->SetFakeRatio(1.);
2864 CookLabel(track,0.); //For comparison only
2866 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
2867 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
2868 if (track->GetNumberOfClusters()<maxn) continue;
2869 maxn = track->GetNumberOfClusters();
2876 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
2877 delete array->RemoveAt(itrack);
2881 if (!besttrack) return;
2884 //take errors of best track as a reference
2885 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
2886 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
2887 for (Int_t i=0;i<6;i++) {
2888 if (besttrack->GetClIndex(i)>0){
2889 erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2890 errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2891 ny[i] = besttrack->GetNy(i);
2892 nz[i] = besttrack->GetNz(i);
2896 // calculate normalized chi2
2898 Float_t * chi2 = new Float_t[entries];
2899 Int_t * index = new Int_t[entries];
2900 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2901 for (Int_t itrack=0;itrack<entries;itrack++){
2902 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2904 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
2905 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
2906 chi2[itrack] = track->GetChi2MIP(0);
2908 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
2909 delete array->RemoveAt(itrack);
2915 TMath::Sort(entries,chi2,index,kFALSE);
2916 besttrack = (AliITStrackMI*)array->At(index[0]);
2917 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
2918 for (Int_t i=0;i<6;i++){
2919 if (besttrack->GetClIndex(i)>0){
2920 erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2921 errz[i] = besttrack->GetSigmaZ(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2922 ny[i] = besttrack->GetNy(i);
2923 nz[i] = besttrack->GetNz(i);
2928 // calculate one more time with updated normalized errors
2929 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2930 for (Int_t itrack=0;itrack<entries;itrack++){
2931 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2933 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2934 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
2935 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
2938 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
2939 delete array->RemoveAt(itrack);
2944 entries = array->GetEntriesFast();
2948 TObjArray * newarray = new TObjArray();
2949 TMath::Sort(entries,chi2,index,kFALSE);
2950 besttrack = (AliITStrackMI*)array->At(index[0]);
2953 for (Int_t i=0;i<6;i++){
2954 if (besttrack->GetNz(i)>0&&besttrack->GetNy(i)>0){
2955 erry[i] = besttrack->GetSigmaY(i); erry[i+6] = besttrack->GetSigmaY(i+6);
2956 errz[i] = besttrack->GetSigmaZ(i); errz[i+6] = besttrack->GetSigmaZ(i+6);
2957 ny[i] = besttrack->GetNy(i);
2958 nz[i] = besttrack->GetNz(i);
2961 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
2962 Float_t minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
2963 Float_t minn = besttrack->GetNumberOfClusters()-3;
2965 for (Int_t i=0;i<entries;i++){
2966 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
2967 if (!track) continue;
2968 if (accepted>maxcut) break;
2969 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
2970 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
2971 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
2972 delete array->RemoveAt(index[i]);
2976 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
2977 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
2978 if (!shortbest) accepted++;
2980 newarray->AddLast(array->RemoveAt(index[i]));
2981 for (Int_t i=0;i<6;i++){
2983 erry[i] = track->GetSigmaY(i); erry[i+6] = track->GetSigmaY(i+6);
2984 errz[i] = track->GetSigmaZ(i); errz[i] = track->GetSigmaZ(i+6);
2985 ny[i] = track->GetNy(i);
2986 nz[i] = track->GetNz(i);
2991 delete array->RemoveAt(index[i]);
2995 delete fTrackHypothesys.RemoveAt(esdindex);
2996 fTrackHypothesys.AddAt(newarray,esdindex);
3000 delete fTrackHypothesys.RemoveAt(esdindex);
3006 //------------------------------------------------------------------------
3007 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3009 //-------------------------------------------------------------
3010 // try to find best hypothesy
3011 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3012 //-------------------------------------------------------------
3013 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3014 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3015 if (!array) return 0;
3016 Int_t entries = array->GetEntriesFast();
3017 if (!entries) return 0;
3018 Float_t minchi2 = 100000;
3019 AliITStrackMI * besttrack=0;
3021 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3022 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3023 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3024 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3026 for (Int_t i=0;i<entries;i++){
3027 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3028 if (!track) continue;
3029 Float_t sigmarfi,sigmaz;
3030 GetDCASigma(track,sigmarfi,sigmaz);
3031 track->SetDnorm(0,sigmarfi);
3032 track->SetDnorm(1,sigmaz);
3034 track->SetChi2MIP(1,1000000);
3035 track->SetChi2MIP(2,1000000);
3036 track->SetChi2MIP(3,1000000);
3039 backtrack = new(backtrack) AliITStrackMI(*track);
3040 if (track->GetConstrain()) {
3041 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3042 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3043 backtrack->ResetCovariance(10.);
3045 backtrack->ResetCovariance(10.);
3047 backtrack->ResetClusters();
3049 Double_t x = original->GetX();
3050 if (!RefitAt(x,backtrack,track)) continue;
3052 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3053 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3054 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3055 track->SetChi22(GetMatchingChi2(backtrack,original));
3057 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3058 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3059 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3062 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3064 //forward track - without constraint
3065 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3066 forwardtrack->ResetClusters();
3068 RefitAt(x,forwardtrack,track);
3069 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3070 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3071 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3073 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3074 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3075 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3076 forwardtrack->SetD(0,track->GetD(0));
3077 forwardtrack->SetD(1,track->GetD(1));
3080 AliITSRecPoint* clist[6];
3081 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3082 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3085 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3086 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3087 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3088 track->SetChi2MIP(3,1000);
3091 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3093 for (Int_t ichi=0;ichi<5;ichi++){
3094 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3096 if (chi2 < minchi2){
3097 //besttrack = new AliITStrackMI(*forwardtrack);
3099 besttrack->SetLabel(track->GetLabel());
3100 besttrack->SetFakeRatio(track->GetFakeRatio());
3102 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3103 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3104 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3108 delete forwardtrack;
3110 for (Int_t i=0;i<entries;i++){
3111 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3112 if (!track) continue;
3114 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3115 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3116 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3117 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3118 delete array->RemoveAt(i);
3128 SortTrackHypothesys(esdindex,checkmax,1);
3129 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3130 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3131 besttrack = (AliITStrackMI*)array->At(0);
3132 if (!besttrack) return 0;
3133 besttrack->SetChi2MIP(8,0);
3134 fBestTrackIndex[esdindex]=0;
3135 entries = array->GetEntriesFast();
3136 AliITStrackMI *longtrack =0;
3138 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3139 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3140 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3141 if (!track->GetConstrain()) continue;
3142 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3143 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3144 if (track->GetChi2MIP(0)>4.) continue;
3145 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3148 //if (longtrack) besttrack=longtrack;
3151 AliITSRecPoint * clist[6];
3152 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3153 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3154 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3155 RegisterClusterTracks(besttrack,esdindex);
3162 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3163 if (sharedtrack>=0){
3165 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3167 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3173 if (shared>2.5) return 0;
3174 if (shared>1.0) return besttrack;
3176 // Don't sign clusters if not gold track
3178 if (!besttrack->IsGoldPrimary()) return besttrack;
3179 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3181 if (fConstraint[fPass]){
3185 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3186 for (Int_t i=0;i<6;i++){
3187 Int_t index = besttrack->GetClIndex(i);
3188 if (index<=0) continue;
3189 Int_t ilayer = (index & 0xf0000000) >> 28;
3190 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3191 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3193 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3194 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3195 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3196 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3197 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3198 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3200 Bool_t cansign = kTRUE;
3201 for (Int_t itrack=0;itrack<entries; itrack++){
3202 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3203 if (!track) continue;
3204 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3205 if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3211 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3212 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3213 if (!c->IsUsed()) c->Use();
3219 //------------------------------------------------------------------------
3220 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3223 // get "best" hypothesys
3226 Int_t nentries = itsTracks.GetEntriesFast();
3227 for (Int_t i=0;i<nentries;i++){
3228 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3229 if (!track) continue;
3230 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3231 if (!array) continue;
3232 if (array->GetEntriesFast()<=0) continue;
3234 AliITStrackMI* longtrack=0;
3236 Float_t maxchi2=1000;
3237 for (Int_t j=0;j<array->GetEntriesFast();j++){
3238 AliITStrackMI* track = (AliITStrackMI*)array->At(j);
3239 if (!track) continue;
3240 if (track->GetGoldV0()) {
3241 longtrack = track; //gold V0 track taken
3244 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3245 Float_t chi2 = track->GetChi2MIP(0);
3247 if (!track->GetGoldV0()&&track->GetConstrain()==kFALSE) chi2+=5;
3249 if (track->GetNumberOfClusters()+track->GetNDeadZone()>minn) maxchi2 = track->GetChi2MIP(0);
3251 if (chi2 > maxchi2) continue;
3252 minn= track->GetNumberOfClusters()+track->GetNDeadZone();
3259 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3260 if (!longtrack) {longtrack = besttrack;}
3261 else besttrack= longtrack;
3265 AliITSRecPoint * clist[6];
3266 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3268 track->SetNUsed(shared);
3269 track->SetNSkipped(besttrack->GetNSkipped());
3270 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3272 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3276 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3277 //if (sharedtrack==-1) sharedtrack=0;
3278 if (sharedtrack>=0) {
3279 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3282 if (besttrack&&fAfterV0) {
3283 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3285 if (besttrack&&fConstraint[fPass])
3286 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3287 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3288 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3289 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3295 //------------------------------------------------------------------------
3296 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3297 //--------------------------------------------------------------------
3298 //This function "cooks" a track label. If label<0, this track is fake.
3299 //--------------------------------------------------------------------
3302 if ( track->GetESDtrack()) tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3304 track->SetChi2MIP(9,0);
3306 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3307 Int_t cindex = track->GetClusterIndex(i);
3308 Int_t l=(cindex & 0xf0000000) >> 28;
3309 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3311 for (Int_t ind=0;ind<3;ind++){
3313 if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3315 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3318 Int_t nclusters = track->GetNumberOfClusters();
3319 if (nclusters > 0) //PH Some tracks don't have any cluster
3320 track->SetFakeRatio(double(nwrong)/double(nclusters));
3322 if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3324 track->SetLabel(tpcLabel);
3328 //------------------------------------------------------------------------
3329 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3334 //AliITSRecPoint * clist[6];
3335 // Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3338 track->SetChi2MIP(9,0);
3339 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3340 Int_t cindex = track->GetClusterIndex(i);
3341 Int_t l=(cindex & 0xf0000000) >> 28;
3342 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3343 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3345 for (Int_t ind=0;ind<3;ind++){
3346 if (cl->GetLabel(ind)==lab) isWrong=0;
3348 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3350 //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue; //shared track
3351 //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3352 //if (l<4&& !(cl->GetType()==1)) continue;
3353 dedx[accepted]= track->GetSampledEdx(i);
3354 //dedx[accepted]= track->fNormQ[l];
3362 TMath::Sort(accepted,dedx,indexes,kFALSE);
3365 Double_t nl=low*accepted, nu =up*accepted;
3367 Float_t sumweight =0;
3368 for (Int_t i=0; i<accepted; i++) {
3370 if (i<nl+0.1) weight = TMath::Max(1.-(nl-i),0.);
3371 if (i>nu-1) weight = TMath::Max(nu-i,0.);
3372 sumamp+= dedx[indexes[i]]*weight;
3375 track->SetdEdx(sumamp/sumweight);
3377 //------------------------------------------------------------------------
3378 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3381 if (fCoefficients) delete []fCoefficients;
3382 fCoefficients = new Float_t[ntracks*48];
3383 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3385 //------------------------------------------------------------------------
3386 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3392 Float_t theta = track->GetTgl();
3393 Float_t phi = track->GetSnp();
3394 phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
3395 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3396 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3398 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3399 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3401 chi2+=0.5*TMath::Min(delta/2,2.);
3402 chi2+=2.*cluster->GetDeltaProbability();
3405 track->SetNy(layer,ny);
3406 track->SetNz(layer,nz);
3407 track->SetSigmaY(layer,erry);
3408 track->SetSigmaZ(layer, errz);
3409 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3410 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/(1.- track->GetSnp()*track->GetSnp())));
3414 //------------------------------------------------------------------------
3415 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3420 Int_t layer = (index & 0xf0000000) >> 28;
3421 track->SetClIndex(layer, index);
3422 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3423 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3424 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3425 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3429 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
3431 //Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);printf("gtr %f %f %f\n",trxyz[0],trxyz[1],trxyz[2]);printf("gcl %f %f %f\n",clxyz[0],clxyz[1],clxyz[2]);
3434 // Take into account the mis-alignment
3435 Double_t x=track->GetX()+cl->GetX();
3436 if (!track->PropagateTo(x,0.,0.)) return 0;
3439 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3440 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3442 return track->UpdateMI(&c,chi2,index);
3445 //------------------------------------------------------------------------
3446 void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3449 //DCA sigmas parameterization
3450 //to be paramterized using external parameters in future
3453 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3454 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3456 //------------------------------------------------------------------------
3457 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3461 Int_t entries = ClusterArray->GetEntriesFast();
3462 if (entries<4) return;
3463 AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3464 Int_t layer = cluster->GetLayer();
3465 if (layer>1) return;
3467 Int_t ncandidates=0;
3468 Float_t r = (layer>0)? 7:4;
3470 for (Int_t i=0;i<entries;i++){
3471 AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3472 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3473 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3474 index[ncandidates] = i; //candidate to belong to delta electron track
3476 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3477 cl0->SetDeltaProbability(1);
3483 for (Int_t i=0;i<ncandidates;i++){
3484 AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3485 if (cl0->GetDeltaProbability()>0.8) continue;
3488 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3489 sumy=sumz=sumy2=sumyz=sumw=0.0;
3490 for (Int_t j=0;j<ncandidates;j++){
3492 AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3494 Float_t dz = cl0->GetZ()-cl1->GetZ();
3495 Float_t dy = cl0->GetY()-cl1->GetY();
3496 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3497 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3498 y[ncl] = cl1->GetY();
3499 z[ncl] = cl1->GetZ();
3500 sumy+= y[ncl]*weight;
3501 sumz+= z[ncl]*weight;
3502 sumy2+=y[ncl]*y[ncl]*weight;
3503 sumyz+=y[ncl]*z[ncl]*weight;
3508 if (ncl<4) continue;
3509 Float_t det = sumw*sumy2 - sumy*sumy;
3511 if (TMath::Abs(det)>0.01){
3512 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3513 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3514 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3517 Float_t z0 = sumyz/sumy;
3518 delta = TMath::Abs(cl0->GetZ()-z0);
3521 cl0->SetDeltaProbability(1-20.*delta);
3525 //------------------------------------------------------------------------
3526 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3530 track->UpdateESDtrack(flags);
3531 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3532 if (oldtrack) delete oldtrack;
3533 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3534 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3535 printf("Problem\n");
3538 //------------------------------------------------------------------------
3539 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3541 // Get nearest upper layer close to the point xr.
3542 // rough approximation
3544 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3545 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3547 for (Int_t i=0;i<6;i++){
3548 if (radius<kRadiuses[i]){
3555 //------------------------------------------------------------------------
3556 void AliITStrackerMI::UpdateTPCV0(AliESDEvent *event){
3558 //try to update, or reject TPC V0s
3560 Int_t nv0s = event->GetNumberOfV0s();
3561 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3563 for (Int_t i=0;i<nv0s;i++){
3564 AliESDv0 * vertex = event->GetV0(i);
3565 Int_t ip = vertex->GetIndex(0);
3566 Int_t im = vertex->GetIndex(1);
3568 TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3569 TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3570 AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3571 AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3575 if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3576 if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3577 if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3582 if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3583 if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3584 if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3587 if (vertex->GetStatus()==-100) continue;
3589 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
3590 Int_t clayer = GetNearestLayer(xrp); //I.B.
3591 vertex->SetNBefore(clayer); //
3592 vertex->SetChi2Before(9*clayer); //
3593 vertex->SetNAfter(6-clayer); //
3594 vertex->SetChi2After(0); //
3596 if (clayer >1 ){ // calculate chi2 before vertex
3597 Float_t chi2p = 0, chi2m=0;
3600 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3601 if (trackp->GetClIndex(ilayer)>0){
3602 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3603 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3614 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3615 if (trackm->GetClIndex(ilayer)>0){
3616 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3617 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3626 vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3627 if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex
3630 if (clayer < 5 ){ // calculate chi2 after vertex
3631 Float_t chi2p = 0, chi2m=0;
3633 if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3634 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3635 if (trackp->GetClIndex(ilayer)>0){
3636 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3637 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3647 if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3648 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3649 if (trackm->GetClIndex(ilayer)>0){
3650 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3651 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3660 vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3661 if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS
3666 //------------------------------------------------------------------------
3667 void AliITStrackerMI::FindV02(AliESDEvent *event)
3672 // Cuts on DCA - R dependent
3673 // max distance DCA between 2 tracks cut
3674 // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3676 const Float_t kMaxDist0 = 0.1;
3677 const Float_t kMaxDist1 = 0.1;
3678 const Float_t kMaxDist = 1;
3679 const Float_t kMinPointAngle = 0.85;
3680 const Float_t kMinPointAngle2 = 0.99;
3681 const Float_t kMinR = 0.5;
3682 const Float_t kMaxR = 220;
3683 //const Float_t kCausality0Cut = 0.19;
3684 //const Float_t kLikelihood01Cut = 0.25;
3685 //const Float_t kPointAngleCut = 0.9996;
3686 const Float_t kCausality0Cut = 0.19;
3687 const Float_t kLikelihood01Cut = 0.45;
3688 const Float_t kLikelihood1Cut = 0.5;
3689 const Float_t kCombinedCut = 0.55;
3693 TTreeSRedirector &cstream = *fDebugStreamer;
3694 Int_t ntracks = event->GetNumberOfTracks();
3695 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3696 fOriginal.Expand(ntracks);
3697 fTrackHypothesys.Expand(ntracks);
3698 fBestHypothesys.Expand(ntracks);
3700 AliHelix * helixes = new AliHelix[ntracks+2];
3701 TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain
3702 TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain
3703 TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain
3704 Bool_t * forbidden = new Bool_t [ntracks+2];
3705 Int_t *itsmap = new Int_t [ntracks+2];
3706 Float_t *dist = new Float_t[ntracks+2];
3707 Float_t *normdist0 = new Float_t[ntracks+2];
3708 Float_t *normdist1 = new Float_t[ntracks+2];
3709 Float_t *normdist = new Float_t[ntracks+2];
3710 Float_t *norm = new Float_t[ntracks+2];
3711 Float_t *maxr = new Float_t[ntracks+2];
3712 Float_t *minr = new Float_t[ntracks+2];
3713 Float_t *minPointAngle= new Float_t[ntracks+2];
3715 AliV0 *pvertex = new AliV0;
3716 AliITStrackMI * dummy= new AliITStrackMI;
3718 AliITStrackMI trackat0; //temporary track for DCA calculation
3720 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3722 // make ITS - ESD map
3724 for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3725 itsmap[itrack] = -1;
3726 forbidden[itrack] = kFALSE;
3727 maxr[itrack] = kMaxR;
3728 minr[itrack] = kMinR;
3729 minPointAngle[itrack] = kMinPointAngle;
3731 for (Int_t itrack=0;itrack<nitstracks;itrack++){
3732 AliITStrackMI * original = (AliITStrackMI*)(fOriginal.At(itrack));
3733 Int_t esdindex = original->GetESDtrack()->GetID();
3734 itsmap[esdindex] = itrack;
3737 // create ITS tracks from ESD tracks if not done before
3739 for (Int_t itrack=0;itrack<ntracks;itrack++){
3740 if (itsmap[itrack]>=0) continue;
3741 AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3742 //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3743 //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3744 tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP()); //I.B.
3745 if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3746 // tracks which can reach inner part of ITS
3747 // propagate track to outer its volume - with correction for material
3748 CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
3750 itsmap[itrack] = nitstracks;
3751 fOriginal.AddAt(tpctrack,nitstracks);
3755 // fill temporary arrays
3757 for (Int_t itrack=0;itrack<ntracks;itrack++){
3758 AliESDtrack * esdtrack = event->GetTrack(itrack);
3759 Int_t itsindex = itsmap[itrack];
3760 AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
3761 if (!original) continue;
3762 AliITStrackMI *bestConst = 0;
3763 AliITStrackMI *bestLong = 0;
3764 AliITStrackMI *best = 0;
3767 TObjArray * array = (TObjArray*) fTrackHypothesys.At(itsindex);
3768 Int_t hentries = (array==0) ? 0 : array->GetEntriesFast();
3769 // Get best track with vertex constrain
3770 for (Int_t ih=0;ih<hentries;ih++){
3771 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3772 if (!trackh->GetConstrain()) continue;
3773 if (!bestConst) bestConst = trackh;
3774 if (trackh->GetNumberOfClusters()>5.0){
3775 bestConst = trackh; // full track - with minimal chi2
3778 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()) continue;
3782 // Get best long track without vertex constrain and best track without vertex constrain
3783 for (Int_t ih=0;ih<hentries;ih++){
3784 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3785 if (trackh->GetConstrain()) continue;
3786 if (!best) best = trackh;
3787 if (!bestLong) bestLong = trackh;
3788 if (trackh->GetNumberOfClusters()>5.0){
3789 bestLong = trackh; // full track - with minimal chi2
3792 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()) continue;
3797 bestLong = original;
3799 //I.B. trackat0 = *bestLong;
3800 new (&trackat0) AliITStrackMI(*bestLong);
3801 Double_t xx,yy,zz,alpha;
3802 bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
3803 alpha = TMath::ATan2(yy,xx);
3804 trackat0.Propagate(alpha,0);
3805 // calculate normalized distances to the vertex
3807 Float_t ptfac = (1.+100.*TMath::Abs(trackat0.GetC()));
3808 if ( bestLong->GetNumberOfClusters()>3 ){
3809 dist[itsindex] = trackat0.GetY();
3810 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3811 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3812 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3813 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3815 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
3816 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
3817 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
3819 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
3820 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
3824 if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
3825 dist[itsindex] = bestConst->GetD(0);
3826 norm[itsindex] = bestConst->GetDnorm(0);
3827 normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3828 normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
3829 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3831 dist[itsindex] = trackat0.GetY();
3832 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
3833 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
3834 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
3835 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
3836 if (TMath::Abs(trackat0.GetTgl())>1.05){
3837 if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
3838 if (normdist[itsindex]>3) {
3839 minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
3845 //-----------------------------------------------------------
3846 //Forbid primary track candidates -
3848 //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
3849 //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
3850 //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
3851 //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
3852 //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
3853 //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
3854 //-----------------------------------------------------------
3856 if (bestLong->GetNumberOfClusters()<4 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5) forbidden[itsindex]=kTRUE;
3857 if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5) forbidden[itsindex]=kTRUE;
3858 if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
3859 if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0) forbidden[itsindex]=kTRUE;
3860 if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2) forbidden[itsindex]=kTRUE;
3861 if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1) forbidden[itsindex]=kTRUE;
3862 if (bestConst->GetNormChi2(0)<2.5) {
3863 minPointAngle[itsindex]= 0.9999;
3864 maxr[itsindex] = 10;
3868 //forbid daughter kink candidates
3870 if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
3871 Bool_t isElectron = kTRUE;
3872 Bool_t isProton = kTRUE;
3874 esdtrack->GetESDpid(pid);
3875 for (Int_t i=1;i<5;i++){
3876 if (pid[0]<pid[i]) isElectron= kFALSE;
3877 if (pid[4]<pid[i]) isProton= kFALSE;
3880 forbidden[itsindex]=kFALSE;
3881 normdist[itsindex]*=-1;
3884 if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
3885 normdist[itsindex]*=-1;
3889 // Causality cuts in TPC volume
3891 if (esdtrack->GetTPCdensity(0,10) >0.6) maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
3892 if (esdtrack->GetTPCdensity(10,30)>0.6) maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
3893 if (esdtrack->GetTPCdensity(20,40)>0.6) maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
3894 if (esdtrack->GetTPCdensity(30,50)>0.6) maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
3896 if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
3902 "Tr1.="<<((bestConst)? bestConst:dummy)<<
3904 "Tr3.="<<&trackat0<<
3906 "Dist="<<dist[itsindex]<<
3907 "ND0="<<normdist0[itsindex]<<
3908 "ND1="<<normdist1[itsindex]<<
3909 "ND="<<normdist[itsindex]<<
3910 "Pz="<<primvertex[2]<<
3911 "Forbid="<<forbidden[itsindex]<<
3915 trackarray.AddAt(best,itsindex);
3916 trackarrayc.AddAt(bestConst,itsindex);
3917 trackarrayl.AddAt(bestLong,itsindex);
3918 new (&helixes[itsindex]) AliHelix(*best);
3923 // first iterration of V0 finder
3925 for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
3926 Int_t itrack0 = itsmap[iesd0];
3927 if (forbidden[itrack0]) continue;
3928 AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
3929 if (!btrack0) continue;
3930 if (btrack0->GetSign()>0) continue;
3931 AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
3933 for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
3934 Int_t itrack1 = itsmap[iesd1];
3935 if (forbidden[itrack1]) continue;
3937 AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
3938 if (!btrack1) continue;
3939 if (btrack1->GetSign()<0) continue;
3940 Bool_t isGold = kFALSE;
3941 if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
3944 AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
3945 AliHelix &h1 = helixes[itrack0];
3946 AliHelix &h2 = helixes[itrack1];
3948 // find linear distance
3953 Double_t phase[2][2],radius[2];
3954 Int_t points = h1.GetRPHIintersections(h2, phase, radius);
3955 if (points==0) continue;
3956 Double_t delta[2]={1000000,1000000};
3958 h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
3960 if (radius[1]<rmin) rmin = radius[1];
3961 h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
3963 rmin = TMath::Sqrt(rmin);
3964 Double_t distance = 0;
3965 Double_t radiusC = 0;
3967 if (points==1 || delta[0]<delta[1]){
3968 distance = TMath::Sqrt(delta[0]);
3969 radiusC = TMath::Sqrt(radius[0]);
3971 distance = TMath::Sqrt(delta[1]);
3972 radiusC = TMath::Sqrt(radius[1]);
3975 if (radiusC<TMath::Max(minr[itrack0],minr[itrack1])) continue;
3976 if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1])) continue;
3977 Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
3978 if (distance>maxDist) continue;
3979 Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
3980 if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
3983 // Double_t distance = TestV0(h1,h2,pvertex,rmin);
3985 // if (distance>maxDist) continue;
3986 // if (pvertex->GetRr()<kMinR) continue;
3987 // if (pvertex->GetRr()>kMaxR) continue;
3988 AliITStrackMI * track0=btrack0;
3989 AliITStrackMI * track1=btrack1;
3990 // if (pvertex->GetRr()<3.5){
3992 //use longest tracks inside the pipe
3993 track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
3994 track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
3998 pvertex->SetParamN(*track0);
3999 pvertex->SetParamP(*track1);
4000 pvertex->Update(primvertex);
4001 pvertex->SetClusters(track0->ClIndex(),track1->ClIndex()); // register clusters
4003 if (pvertex->GetRr()<kMinR) continue;
4004 if (pvertex->GetRr()>kMaxR) continue;
4005 if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4006 //Bo: if (pvertex->GetDist2()>maxDist) continue;
4007 if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4008 //Bo: pvertex->SetLab(0,track0->GetLabel());
4009 //Bo: pvertex->SetLab(1,track1->GetLabel());
4010 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4011 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4013 AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4014 AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4018 TObjArray * array0b = (TObjArray*)fBestHypothesys.At(itrack0);
4019 if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4020 fCurrentEsdTrack = itrack0;
4021 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4023 TObjArray * array1b = (TObjArray*)fBestHypothesys.At(itrack1);
4024 if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4025 fCurrentEsdTrack = itrack1;
4026 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4029 AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4030 AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4031 AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4032 AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4034 Float_t minchi2before0=16;
4035 Float_t minchi2before1=16;
4036 Float_t minchi2after0 =16;
4037 Float_t minchi2after1 =16;
4038 Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
4039 Int_t maxLayer = GetNearestLayer(xrp); //I.B.
4041 if (array0b) for (Int_t i=0;i<5;i++){
4042 // best track after vertex
4043 AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4044 if (!btrack) continue;
4045 if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4046 // if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4047 if (btrack->GetX()<pvertex->GetRr()-2.) {
4048 if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4051 if (maxLayer<3){ // take prim vertex as additional measurement
4052 if (normdist[itrack0]>0 && htrackc0){
4053 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4055 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4059 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4061 if (!btrack->GetClIndex(ilayer)){
4065 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4066 for (Int_t itrack=0;itrack<4;itrack++){
4067 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4068 sumchi2+=18.; //shared cluster
4072 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4073 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4077 if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4079 continue; //safety space - Geo manager will give exact layer
4082 minchi2after0 = btrack->GetNormChi2(i);
4085 if (array1b) for (Int_t i=0;i<5;i++){
4086 // best track after vertex
4087 AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4088 if (!btrack) continue;
4089 if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4090 // if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4091 if (btrack->GetX()<pvertex->GetRr()-2){
4092 if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4095 if (maxLayer<3){ // take prim vertex as additional measurement
4096 if (normdist[itrack1]>0 && htrackc1){
4097 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4099 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4103 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4105 if (!btrack->GetClIndex(ilayer)){
4109 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4110 for (Int_t itrack=0;itrack<4;itrack++){
4111 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4112 sumchi2+=18.; //shared cluster
4116 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4117 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4121 if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4123 continue; //safety space - Geo manager will give exact layer
4126 minchi2after1 = btrack->GetNormChi2(i);
4130 // position resolution - used for DCA cut
4131 Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4132 (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4133 (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4134 sigmad =TMath::Sqrt(sigmad)+0.04;
4135 if (pvertex->GetRr()>50){
4136 Double_t cov0[15],cov1[15];
4137 track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4138 track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4139 sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4140 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4141 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4142 sigmad =TMath::Sqrt(sigmad)+0.05;
4146 vertex2.SetParamN(*track0b);
4147 vertex2.SetParamP(*track1b);
4148 vertex2.Update(primvertex);
4149 //Bo: if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4150 if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4151 pvertex->SetParamN(*track0b);
4152 pvertex->SetParamP(*track1b);
4153 pvertex->Update(primvertex);
4154 pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex()); // register clusters
4155 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4156 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4158 pvertex->SetDistSigma(sigmad);
4159 //Bo: pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4160 pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4162 // define likelihhod and causalities
4164 Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4166 Float_t fnorm0 = normdist[itrack0];
4167 if (fnorm0<0) fnorm0*=-3;
4168 Float_t fnorm1 = normdist[itrack1];
4169 if (fnorm1<0) fnorm1*=-3;
4170 if (pvertex->GetAnglep()[2]>0.1 || (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 || pvertex->GetRr()<3){
4171 pb0 = TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4172 pb1 = TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4174 pvertex->SetChi2Before(normdist[itrack0]);
4175 pvertex->SetChi2After(normdist[itrack1]);
4176 pvertex->SetNAfter(0);
4177 pvertex->SetNBefore(0);
4179 pvertex->SetChi2Before(minchi2before0);
4180 pvertex->SetChi2After(minchi2before1);
4181 if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4182 pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4183 pb1 = TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4185 pvertex->SetNAfter(maxLayer);
4186 pvertex->SetNBefore(maxLayer);
4188 if (pvertex->GetRr()<90){
4189 pa0 *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4190 pa1 *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4192 if (pvertex->GetRr()<20){
4193 pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4194 pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4197 pvertex->SetCausality(pb0,pb1,pa0,pa1);
4199 // Likelihood calculations - apply cuts
4201 Bool_t v0OK = kTRUE;
4202 Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4203 p12 += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4204 p12 = TMath::Sqrt(p12); // "mean" momenta
4205 Float_t sigmap0 = 0.0001+0.001/(0.1+pvertex->GetRr());
4206 Float_t sigmap = 0.5*sigmap0*(0.6+0.4*p12); // "resolution: of point angle - as a function of radius and momenta
4208 Float_t causalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4209 Float_t causalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4210 TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4212 //Bo: Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4213 Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4214 Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4216 Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4217 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4218 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4219 0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4221 if (causalityA<kCausality0Cut) v0OK = kFALSE;
4222 if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut) v0OK = kFALSE;
4223 if (likelihood1<kLikelihood1Cut) v0OK = kFALSE;
4224 if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4229 Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4231 "Tr0.="<<track0<< //best without constrain
4232 "Tr1.="<<track1<< //best without constrain
4233 "Tr0B.="<<track0b<< //best without constrain after vertex
4234 "Tr1B.="<<track1b<< //best without constrain after vertex
4235 "Tr0C.="<<htrackc0<< //best with constrain if exist
4236 "Tr1C.="<<htrackc1<< //best with constrain if exist
4237 "Tr0L.="<<track0l<< //longest best
4238 "Tr1L.="<<track1l<< //longest best
4239 "Esd0.="<<track0->GetESDtrack()<< // esd track0 params
4240 "Esd1.="<<track1->GetESDtrack()<< // esd track1 params
4241 "V0.="<<pvertex<< //vertex properties
4242 "V0b.="<<&vertex2<< //vertex properties at "best" track
4243 "ND0="<<normdist[itrack0]<< //normalize distance for track0
4244 "ND1="<<normdist[itrack1]<< //normalize distance for track1
4246 // "RejectBase="<<rejectBase<< //rejection in First itteration
4252 //if (rejectBase) continue;
4254 pvertex->SetStatus(0);
4255 // if (rejectBase) {
4256 // pvertex->SetStatus(-100);
4258 if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4259 //Bo: pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4260 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4261 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4263 // AliV0vertex vertexjuri(*track0,*track1);
4264 // vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4265 // event->AddV0(&vertexjuri);
4266 pvertex->SetStatus(100);
4268 pvertex->SetOnFlyStatus(kTRUE);
4269 pvertex->ChangeMassHypothesis(kK0Short);
4270 event->AddV0(pvertex);
4276 // delete temporary arrays
4279 delete[] minPointAngle;
4291 //------------------------------------------------------------------------
4292 void AliITStrackerMI::RefitV02(AliESDEvent *event)
4295 //try to refit V0s in the third path of the reconstruction
4297 TTreeSRedirector &cstream = *fDebugStreamer;
4299 Int_t nv0s = event->GetNumberOfV0s();
4300 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4302 for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4303 AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4304 if (!v0mi) continue;
4305 Int_t itrack0 = v0mi->GetIndex(0);
4306 Int_t itrack1 = v0mi->GetIndex(1);
4307 AliESDtrack *esd0 = event->GetTrack(itrack0);
4308 AliESDtrack *esd1 = event->GetTrack(itrack1);
4309 if (!esd0||!esd1) continue;
4310 AliITStrackMI tpc0(*esd0);
4311 AliITStrackMI tpc1(*esd1);
4312 Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4313 Double_t alpha =TMath::ATan2(y,x); //I.B.
4314 if (v0mi->GetRr()>85){
4315 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4316 v0temp.SetParamN(tpc0);
4317 v0temp.SetParamP(tpc1);
4318 v0temp.Update(primvertex);
4319 if (kFALSE) cstream<<"Refit"<<
4321 "V0refit.="<<&v0temp<<
4325 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4326 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4327 v0mi->SetParamN(tpc0);
4328 v0mi->SetParamP(tpc1);
4329 v0mi->Update(primvertex);
4334 if (v0mi->GetRr()>35){
4335 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4336 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4337 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4338 v0temp.SetParamN(tpc0);
4339 v0temp.SetParamP(tpc1);
4340 v0temp.Update(primvertex);
4341 if (kFALSE) cstream<<"Refit"<<
4343 "V0refit.="<<&v0temp<<
4347 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4348 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4349 v0mi->SetParamN(tpc0);
4350 v0mi->SetParamP(tpc1);
4351 v0mi->Update(primvertex);
4356 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4357 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4358 // if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4359 if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4360 v0temp.SetParamN(tpc0);
4361 v0temp.SetParamP(tpc1);
4362 v0temp.Update(primvertex);
4363 if (kFALSE) cstream<<"Refit"<<
4365 "V0refit.="<<&v0temp<<
4369 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4370 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4371 v0mi->SetParamN(tpc0);
4372 v0mi->SetParamP(tpc1);
4373 v0mi->Update(primvertex);
4378 //------------------------------------------------------------------------
4379 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4380 //--------------------------------------------------------------------
4381 // Fill a look-up table with mean material
4382 //--------------------------------------------------------------------
4386 Double_t point1[3],point2[3];
4387 Double_t phi,cosphi,sinphi,z;
4388 // 0-5 layers, 6 pipe, 7-8 shields
4389 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 7.5,25.0};
4390 Double_t rmax[9]={ 5.5, 7.3,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4392 Int_t ifirst=0,ilast=0;
4393 if(material.Contains("Pipe")) {
4395 } else if(material.Contains("Shields")) {
4397 } else if(material.Contains("Layers")) {
4400 Error("BuildMaterialLUT","Wrong layer name\n");
4403 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4404 Double_t param[5]={0.,0.,0.,0.,0.};
4405 for (Int_t i=0; i<n; i++) {
4406 phi = 2.*TMath::Pi()*gRandom->Rndm();
4407 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4408 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4409 point1[0] = rmin[imat]*cosphi;
4410 point1[1] = rmin[imat]*sinphi;
4412 point2[0] = rmax[imat]*cosphi;
4413 point2[1] = rmax[imat]*sinphi;
4415 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4416 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4418 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4420 fxOverX0Layer[imat] = param[1];
4421 fxTimesRhoLayer[imat] = param[0]*param[4];
4422 } else if(imat==6) {
4423 fxOverX0Pipe = param[1];
4424 fxTimesRhoPipe = param[0]*param[4];
4425 } else if(imat==7) {
4426 fxOverX0Shield[0] = param[1];
4427 fxTimesRhoShield[0] = param[0]*param[4];
4428 } else if(imat==8) {
4429 fxOverX0Shield[1] = param[1];
4430 fxTimesRhoShield[1] = param[0]*param[4];
4434 printf("%s\n",material.Data());
4435 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4436 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4437 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4438 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4439 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4440 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4441 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4442 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4443 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4447 //------------------------------------------------------------------------
4448 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4449 TString direction) {
4450 //-------------------------------------------------------------------
4451 // Propagate beyond beam pipe and correct for material
4452 // (material budget in different ways according to fUseTGeo value)
4453 //-------------------------------------------------------------------
4455 // Define budget mode:
4456 // 0: material from AliITSRecoParam (hard coded)
4457 // 1: material from TGeo (on the fly)
4458 // 2: material from lut
4459 // 3: material from TGeo (same for all hypotheses)
4472 if(fTrackingPhase.Contains("Clusters2Tracks"))
4473 { mode=3; } else { mode=1; }
4476 if(fTrackingPhase.Contains("Clusters2Tracks"))
4477 { mode=3; } else { mode=2; }
4483 if(fTrackingPhase.Contains("Default")) mode=0;
4485 Int_t index=fCurrentEsdTrack;
4487 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4488 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4489 Double_t xToGo; t->GetLocalXat(rToGo,xToGo);
4491 Double_t xOverX0,x0,lengthTimesMeanDensity;
4492 Bool_t anglecorr=kTRUE;
4496 xOverX0 = AliITSRecoParam::GetdPipe();
4497 x0 = AliITSRecoParam::GetX0Be();
4498 lengthTimesMeanDensity = xOverX0*x0;
4501 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4505 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4506 xOverX0 = fxOverX0Pipe;
4507 lengthTimesMeanDensity = fxTimesRhoPipe;
4510 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4511 if(fxOverX0PipeTrks[index]<0) {
4512 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4513 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4514 (1.-t->GetSnp()*t->GetSnp()));
4515 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4516 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4519 xOverX0 = fxOverX0PipeTrks[index];
4520 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4524 lengthTimesMeanDensity *= dir;
4526 if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4527 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4531 //------------------------------------------------------------------------
4532 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4534 TString direction) {
4535 //-------------------------------------------------------------------
4536 // Propagate beyond SPD or SDD shield and correct for material
4537 // (material budget in different ways according to fUseTGeo value)
4538 //-------------------------------------------------------------------
4540 // Define budget mode:
4541 // 0: material from AliITSRecoParam (hard coded)
4542 // 1: material from TGeo (on the fly)
4543 // 2: material from lut
4544 // 3: material from TGeo (same for all hypotheses)
4557 if(fTrackingPhase.Contains("Clusters2Tracks"))
4558 { mode=3; } else { mode=1; }
4561 if(fTrackingPhase.Contains("Clusters2Tracks"))
4562 { mode=3; } else { mode=2; }
4568 if(fTrackingPhase.Contains("Default")) mode=0;
4570 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4572 Int_t shieldindex=0;
4573 if (shield.Contains("SDD")) { // SDDouter
4574 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4576 } else if (shield.Contains("SPD")) { // SPDouter
4577 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4580 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4583 Double_t xToGo; t->GetLocalXat(rToGo,xToGo);
4585 Int_t index=2*fCurrentEsdTrack+shieldindex;
4587 Double_t xOverX0,x0,lengthTimesMeanDensity;
4588 Bool_t anglecorr=kTRUE;
4592 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4593 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4594 lengthTimesMeanDensity = xOverX0*x0;
4597 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4601 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4602 xOverX0 = fxOverX0Shield[shieldindex];
4603 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4606 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4607 if(fxOverX0ShieldTrks[index]<0) {
4608 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4609 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4610 (1.-t->GetSnp()*t->GetSnp()));
4611 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4612 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4615 xOverX0 = fxOverX0ShieldTrks[index];
4616 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4620 lengthTimesMeanDensity *= dir;
4622 if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4623 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4627 //------------------------------------------------------------------------
4628 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4630 Double_t oldGlobXYZ[3],
4631 TString direction) {
4632 //-------------------------------------------------------------------
4633 // Propagate beyond layer and correct for material
4634 // (material budget in different ways according to fUseTGeo value)
4635 //-------------------------------------------------------------------
4637 // Define budget mode:
4638 // 0: material from AliITSRecoParam (hard coded)
4639 // 1: material from TGeo (on the fly)
4640 // 2: material from lut
4641 // 3: material from TGeo (same for all hypotheses)
4654 if(fTrackingPhase.Contains("Clusters2Tracks"))
4655 { mode=3; } else { mode=1; }
4658 if(fTrackingPhase.Contains("Clusters2Tracks"))
4659 { mode=3; } else { mode=2; }
4665 if(fTrackingPhase.Contains("Default")) mode=0;
4667 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4669 Double_t r=fgLayers[layerindex].GetR();
4670 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4672 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4673 Double_t xToGo; t->GetLocalXat(rToGo,xToGo);
4675 Int_t index=6*fCurrentEsdTrack+layerindex;
4677 // Bring the track beyond the material
4678 if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4679 Double_t globXYZ[3];
4682 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4684 Bool_t anglecorr=kTRUE;
4688 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4689 lengthTimesMeanDensity = xOverX0*x0;
4692 AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4693 if(mparam[1]>900000) return 0;
4695 lengthTimesMeanDensity=mparam[0]*mparam[4];
4699 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4700 xOverX0 = fxOverX0Layer[layerindex];
4701 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4704 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4705 if(fxOverX0LayerTrks[index]<0) {
4706 AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4707 if(mparam[1]>900000) return 0;
4708 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4709 (1.-t->GetSnp()*t->GetSnp()));
4710 xOverX0=mparam[1]/angle;
4711 lengthTimesMeanDensity=mparam[0]*mparam[4]/angle;
4712 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0);
4713 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity);
4715 xOverX0 = fxOverX0LayerTrks[index];
4716 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4720 lengthTimesMeanDensity *= dir;
4722 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4726 //------------------------------------------------------------------------
4727 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4728 //-----------------------------------------------------------------
4729 // Initialize LUT for storing material for each prolonged track
4730 //-----------------------------------------------------------------
4731 fxOverX0PipeTrks = new Float_t[ntracks];
4732 fxTimesRhoPipeTrks = new Float_t[ntracks];
4733 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4734 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4735 fxOverX0LayerTrks = new Float_t[ntracks*6];
4736 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4738 for(Int_t i=0; i<ntracks; i++) {
4739 fxOverX0PipeTrks[i] = -1.;
4740 fxTimesRhoPipeTrks[i] = -1.;
4742 for(Int_t j=0; j<ntracks*2; j++) {
4743 fxOverX0ShieldTrks[j] = -1.;
4744 fxTimesRhoShieldTrks[j] = -1.;
4746 for(Int_t k=0; k<ntracks*6; k++) {
4747 fxOverX0LayerTrks[k] = -1.;
4748 fxTimesRhoLayerTrks[k] = -1.;
4755 //------------------------------------------------------------------------
4756 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4757 //-----------------------------------------------------------------
4758 // Delete LUT for storing material for each prolonged track
4759 //-----------------------------------------------------------------
4760 if(fxOverX0PipeTrks) {
4761 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4763 if(fxOverX0ShieldTrks) {
4764 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4767 if(fxOverX0LayerTrks) {
4768 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4770 if(fxTimesRhoPipeTrks) {
4771 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4773 if(fxTimesRhoShieldTrks) {
4774 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4776 if(fxTimesRhoLayerTrks) {
4777 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4781 //------------------------------------------------------------------------
4782 Int_t AliITStrackerMI::CheckSkipLayer(AliITStrackMI *track,
4783 Int_t ilayer,Int_t idet) const {
4784 //-----------------------------------------------------------------
4785 // This method is used to decide whether to allow a prolongation
4786 // without clusters, because we want to skip the layer.
4787 // In this case the return value is > 0:
4788 // return 1: the user requested to skip a layer
4789 // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
4790 //-----------------------------------------------------------------
4792 if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
4794 if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4795 // check if track will cross SPD outer layer
4796 Double_t phiAtSPD2,zAtSPD2;
4797 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4798 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4804 //------------------------------------------------------------------------
4805 Int_t AliITStrackerMI::CheckDeadZone(/*AliITStrackMI *track,*/
4806 Int_t ilayer,Int_t idet,
4807 Double_t zmin,Double_t zmax/*,Double_t ymin,Double_t ymax*/) const {
4808 //-----------------------------------------------------------------
4809 // This method is used to decide whether to allow a prolongation
4810 // without clusters, because there is a dead zone in the road.
4811 // In this case the return value is > 0:
4812 // return 1: dead zone at z=0,+-7cm in SPD
4813 // return 2: dead area from the OCDB // NOT YET IMPLEMENTED
4814 //-----------------------------------------------------------------
4816 // check dead zones at z=0,+-7cm in the SPD
4817 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4818 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4819 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4820 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4821 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4822 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4823 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4824 for (Int_t i=0; i<3; i++)
4825 if (zmin<zmaxdead[i] && zmax>zmindead[i]) return 1;
4828 // check dead zones from OCDB
4829 if (!AliITSReconstructor::GetRecoParam()->GetUseDeadZonesFromOCDB()) return 0;
4831 if(idet<0) return 0;
4833 // look in OCDB (only entire dead modules for the moment)
4834 if (ilayer==0 || ilayer==1) { // SPD
4835 AliCDBEntry* cdbEntry = AliCDBManager::Instance()->Get("ITS/Calib/SPDDead");
4837 Error("CheckDeadZone","Cannot get CDB entry for SPD\n");
4840 TObjArray* spdEntry = (TObjArray*)cdbEntry->GetObject();
4842 Error("CheckDeadZone","Cannot get CDB entry for SPD\n");
4845 if(ilayer==1) idet += AliITSgeomTGeo::GetNLadders(1)*AliITSgeomTGeo::GetNDetectors(1);
4846 //printf("SPD det: %d\n",idet);
4847 AliITSCalibrationSPD *calibSPD = (AliITSCalibrationSPD*)spdEntry->At(idet);
4848 if (calibSPD->IsBad()) return 2;
4849 } else if (ilayer==2 || ilayer==3) { // SDD
4850 AliCDBEntry* cdbEntry = AliCDBManager::Instance()->Get("ITS/Calib/CalibSDD");
4852 Error("CheckDeadZone","Cannot get CDB entry for SDD\n");
4855 TObjArray* sddEntry = (TObjArray*)cdbEntry->GetObject();
4857 Error("CheckDeadZone","Cannot get CDB entry for SDD\n");
4860 if(ilayer==3) idet += AliITSgeomTGeo::GetNLadders(3)*AliITSgeomTGeo::GetNDetectors(3);
4861 //printf("SDD det: %d\n",idet);
4862 AliITSCalibrationSDD *calibSDD = (AliITSCalibrationSDD*)sddEntry->At(idet);
4863 if (calibSDD->IsDead()) return 2;
4864 } else if (ilayer==4 || ilayer==5) { // SSD
4866 Error("CheckDeadZone","Wrong layer number\n");
4867 if(ilayer==5) idet += AliITSgeomTGeo::GetNLadders(5)*AliITSgeomTGeo::GetNDetectors(5);
4873 //------------------------------------------------------------------------
4874 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4875 AliITStrackMI *track,
4876 Float_t &xloc,Float_t &zloc) const {
4877 //-----------------------------------------------------------------
4878 // Gives position of track in local module ref. frame
4879 //-----------------------------------------------------------------
4884 if(idet<0) return kFALSE;
4886 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4888 Int_t lad = Int_t(idet/ndet) + 1;
4890 Int_t det = idet - (lad-1)*ndet + 1;
4892 Double_t xyzGlob[3],xyzLoc[3];
4894 track->GetXYZ(xyzGlob);
4896 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
4898 xloc = (Float_t)xyzLoc[0];
4899 zloc = (Float_t)xyzLoc[2];
4903 //------------------------------------------------------------------------