1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------------------------
18 // Implementation of the ITS tracker class
19 // It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
20 // and fills with them the ESD
21 // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
22 // Current support and development:
23 // Andrea Dainese, andrea.dainese@lnl.infn.it
24 // dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25 // Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26 // Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
27 //-------------------------------------------------------------------------
31 #include <TTreeStream.h>
32 #include <TDatabasePDG.h>
38 #include "AliESDEvent.h"
39 #include "AliESDtrack.h"
40 #include "AliESDVertex.h"
43 #include "AliITSRecPoint.h"
44 #include "AliITSgeomTGeo.h"
45 #include "AliITSReconstructor.h"
46 #include "AliTrackPointArray.h"
47 #include "AliAlignObj.h"
48 #include "AliITSClusterParam.h"
49 #include "AliCDBManager.h"
50 #include "AliCDBEntry.h"
51 #include "AliITSsegmentation.h"
52 #include "AliITSCalibration.h"
53 #include "AliITSCalibrationSPD.h"
54 #include "AliITSCalibrationSDD.h"
55 #include "AliITSCalibrationSSD.h"
56 #include "AliITSPlaneEff.h"
57 #include "AliITSPlaneEffSPD.h"
58 #include "AliITSPlaneEffSDD.h"
59 #include "AliITSPlaneEffSSD.h"
60 #include "AliITStrackerMI.h"
62 ClassImp(AliITStrackerMI)
64 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
66 AliITStrackerMI::AliITStrackerMI():AliTracker(),
76 fLastLayerToTrackTo(0),
79 fTrackingPhase("Default"),
85 fxTimesRhoPipeTrks(0),
86 fxOverX0ShieldTrks(0),
87 fxTimesRhoShieldTrks(0),
89 fxTimesRhoLayerTrks(0),
96 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
97 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
98 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
100 //------------------------------------------------------------------------
101 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
102 fI(AliITSgeomTGeo::GetNLayers()),
111 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
114 fTrackingPhase("Default"),
120 fxTimesRhoPipeTrks(0),
121 fxOverX0ShieldTrks(0),
122 fxTimesRhoShieldTrks(0),
123 fxOverX0LayerTrks(0),
124 fxTimesRhoLayerTrks(0),
126 fITSChannelStatus(0),
129 //--------------------------------------------------------------------
130 //This is the AliITStrackerMI constructor
131 //--------------------------------------------------------------------
133 AliWarning("\"geom\" is actually a dummy argument !");
139 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
140 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
141 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
143 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
144 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
145 Double_t poff=TMath::ATan2(y,x);
147 Double_t r=TMath::Sqrt(x*x + y*y);
149 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
150 r += TMath::Sqrt(x*x + y*y);
151 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
152 r += TMath::Sqrt(x*x + y*y);
153 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
154 r += TMath::Sqrt(x*x + y*y);
157 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
159 for (Int_t j=1; j<nlad+1; j++) {
160 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
161 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
162 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
164 Double_t txyz[3]={0.};
165 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
166 m.LocalToMaster(txyz,xyz);
167 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
168 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
170 if (phi<0) phi+=TMath::TwoPi();
171 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
173 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
174 new(&det) AliITSdetector(r,phi);
175 // compute the real radius (with misalignment)
176 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
178 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
179 mmisal.LocalToMaster(txyz,xyz);
180 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
181 det.SetRmisal(rmisal);
183 } // end loop on detectors
184 } // end loop on ladders
185 } // end loop on layers
188 fI=AliITSgeomTGeo::GetNLayers();
191 fConstraint[0]=1; fConstraint[1]=0;
193 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
194 AliITSReconstructor::GetRecoParam()->GetYVdef(),
195 AliITSReconstructor::GetRecoParam()->GetZVdef()};
196 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
197 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
198 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
199 SetVertex(xyzVtx,ersVtx);
201 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=AliITSRecoParam::GetLayersNotToSkip(i);
202 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
203 for (Int_t i=0;i<100000;i++){
204 fBestTrackIndex[i]=0;
207 // store positions of centre of SPD modules (in z)
209 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
210 fSPDdetzcentre[0] = tr[2];
211 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
212 fSPDdetzcentre[1] = tr[2];
213 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
214 fSPDdetzcentre[2] = tr[2];
215 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
216 fSPDdetzcentre[3] = tr[2];
218 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
219 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
220 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
224 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
225 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
227 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
229 // only for plane efficiency evaluation
230 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff()) {
231 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
232 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane))
233 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
234 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
235 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
236 else fPlaneEff = new AliITSPlaneEffSSD();
237 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
238 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
239 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
242 //------------------------------------------------------------------------
243 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
245 fBestTrack(tracker.fBestTrack),
246 fTrackToFollow(tracker.fTrackToFollow),
247 fTrackHypothesys(tracker.fTrackHypothesys),
248 fBestHypothesys(tracker.fBestHypothesys),
249 fOriginal(tracker.fOriginal),
250 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
251 fPass(tracker.fPass),
252 fAfterV0(tracker.fAfterV0),
253 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
254 fCoefficients(tracker.fCoefficients),
256 fTrackingPhase(tracker.fTrackingPhase),
257 fUseTGeo(tracker.fUseTGeo),
258 fNtracks(tracker.fNtracks),
259 fxOverX0Pipe(tracker.fxOverX0Pipe),
260 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
262 fxTimesRhoPipeTrks(0),
263 fxOverX0ShieldTrks(0),
264 fxTimesRhoShieldTrks(0),
265 fxOverX0LayerTrks(0),
266 fxTimesRhoLayerTrks(0),
267 fDebugStreamer(tracker.fDebugStreamer),
268 fITSChannelStatus(tracker.fITSChannelStatus),
269 fDetTypeRec(tracker.fDetTypeRec),
270 fPlaneEff(tracker.fPlaneEff) {
274 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
277 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
278 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
281 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
282 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
285 //------------------------------------------------------------------------
286 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
287 //Assignment operator
288 this->~AliITStrackerMI();
289 new(this) AliITStrackerMI(tracker);
292 //------------------------------------------------------------------------
293 AliITStrackerMI::~AliITStrackerMI()
298 if (fCoefficients) delete [] fCoefficients;
299 DeleteTrksMaterialLUT();
300 if (fDebugStreamer) {
301 //fDebugStreamer->Close();
302 delete fDebugStreamer;
304 if(fITSChannelStatus) delete fITSChannelStatus;
305 if(fPlaneEff) delete fPlaneEff;
307 //------------------------------------------------------------------------
308 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
309 //--------------------------------------------------------------------
310 //This function set masks of the layers which must be not skipped
311 //--------------------------------------------------------------------
312 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=l[i];
314 //------------------------------------------------------------------------
315 void AliITStrackerMI::ReadBadFromDetTypeRec() {
316 //--------------------------------------------------------------------
317 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
319 //--------------------------------------------------------------------
321 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
323 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
325 if(!fDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
328 if(fITSChannelStatus) delete fITSChannelStatus;
329 fITSChannelStatus = new AliITSChannelStatus(fDetTypeRec);
331 // ITS detectors and chips
332 Int_t i=0,j=0,k=0,ndet=0;
333 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
334 Int_t nBadDetsPerLayer=0;
335 ndet=AliITSgeomTGeo::GetNDetectors(i);
336 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
337 for (k=1; k<ndet+1; k++) {
338 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
339 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fDetTypeRec);
340 if(det.IsBad()) {nBadDetsPerLayer++;}
341 } // end loop on detectors
342 } // end loop on ladders
343 Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
344 } // end loop on layers
348 //------------------------------------------------------------------------
349 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
350 //--------------------------------------------------------------------
351 //This function loads ITS clusters
352 //--------------------------------------------------------------------
353 TBranch *branch=cTree->GetBranch("ITSRecPoints");
355 Error("LoadClusters"," can't get the branch !\n");
359 static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
360 branch->SetAddress(&clusters);
362 Int_t i=0,j=0,ndet=0;
364 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
365 ndet=fgLayers[i].GetNdetectors();
366 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
367 for (; j<jmax; j++) {
368 if (!cTree->GetEvent(j)) continue;
369 Int_t ncl=clusters->GetEntriesFast();
370 SignDeltas(clusters,GetZ());
373 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
374 detector=c->GetDetectorIndex();
376 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
378 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
381 // add dead zone "virtual" cluster in SPD, if there is a cluster within
382 // zwindow cm from the dead zone
383 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
384 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
385 Int_t lab[4] = {0,0,0,detector};
386 Int_t info[3] = {0,0,i};
387 Float_t q = 0.; // this identifies virtual clusters
388 Float_t hit[5] = {xdead,
390 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
391 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
393 Bool_t local = kTRUE;
394 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
395 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
396 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
397 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
398 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
399 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
400 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
401 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
402 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
403 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
404 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
405 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
406 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
407 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
408 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
409 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
410 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
411 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
412 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
414 } // "virtual" clusters in SPD
418 fgLayers[i].ResetRoad(); //road defined by the cluster density
419 fgLayers[i].SortClusters();
426 //------------------------------------------------------------------------
427 void AliITStrackerMI::UnloadClusters() {
428 //--------------------------------------------------------------------
429 //This function unloads ITS clusters
430 //--------------------------------------------------------------------
431 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
433 //------------------------------------------------------------------------
434 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
435 //--------------------------------------------------------------------
436 // Publishes all pointers to clusters known to the tracker into the
437 // passed object array.
438 // The ownership is not transfered - the caller is not expected to delete
440 //--------------------------------------------------------------------
442 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
443 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
444 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
451 //------------------------------------------------------------------------
452 static Int_t CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
453 //--------------------------------------------------------------------
454 // Correction for the material between the TPC and the ITS
455 //--------------------------------------------------------------------
456 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
457 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
458 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
459 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
460 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
461 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
462 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
463 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
465 Error("CorrectForTPCtoITSDeadZoneMaterial","Track is already in the dead zone !");
471 //------------------------------------------------------------------------
472 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
473 //--------------------------------------------------------------------
474 // This functions reconstructs ITS tracks
475 // The clusters must be already loaded !
476 //--------------------------------------------------------------------
479 fTrackingPhase="Clusters2Tracks";
481 TObjArray itsTracks(15000);
483 fEsd = event; // store pointer to the esd
485 // temporary (for cosmics)
486 if(event->GetVertex()) {
487 TString title = event->GetVertex()->GetTitle();
488 if(title.Contains("cosmics")) {
489 Double_t xyz[3]={GetX(),GetY(),GetZ()};
490 Double_t exyz[3]={0.1,0.1,0.1};
496 {/* Read ESD tracks */
497 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
498 Int_t nentr=event->GetNumberOfTracks();
499 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
501 AliESDtrack *esd=event->GetTrack(nentr);
502 // ---- for debugging:
503 //if(TMath::Abs(esd->GetX()-83.65)<0.1) { FILE *f=fopen("tpc.dat","a"); fprintf(f,"%f %f %f %f %f %f\n",(Float_t)event->GetEventNumberInFile(),(Float_t)TMath::Abs(esd->GetLabel()),(Float_t)esd->GetX(),(Float_t)esd->GetY(),(Float_t)esd->GetZ(),(Float_t)esd->Pt()); fclose(f); }
505 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
506 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
507 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
508 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
511 t=new AliITStrackMI(*esd);
512 } catch (const Char_t *msg) {
513 //Warning("Clusters2Tracks",msg);
517 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
518 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
521 // look at the ESD mass hypothesys !
522 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
524 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
526 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
527 //track - can be V0 according to TPC
529 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
533 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
537 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
541 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
546 t->SetReconstructed(kFALSE);
547 itsTracks.AddLast(t);
548 fOriginal.AddLast(t);
550 } /* End Read ESD tracks */
554 Int_t nentr=itsTracks.GetEntriesFast();
555 fTrackHypothesys.Expand(nentr);
556 fBestHypothesys.Expand(nentr);
557 MakeCoefficients(nentr);
558 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
560 // THE TWO TRACKING PASSES
561 for (fPass=0; fPass<2; fPass++) {
562 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
563 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
564 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
565 if (t==0) continue; //this track has been already tracked
566 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
567 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
568 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
569 if (fConstraint[fPass]) {
570 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
571 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
574 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
575 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
577 ResetTrackToFollow(*t);
580 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
583 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
585 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
586 if (!besttrack) continue;
587 besttrack->SetLabel(tpcLabel);
588 // besttrack->CookdEdx();
590 besttrack->SetFakeRatio(1.);
591 CookLabel(besttrack,0.); //For comparison only
592 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
594 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
596 t->SetReconstructed(kTRUE);
598 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
600 GetBestHypothesysMIP(itsTracks);
601 } // end loop on the two tracking passes
603 if(event->GetNumberOfV0s()>0) UpdateTPCV0(event);
604 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) FindV02(event);
609 Int_t entries = fTrackHypothesys.GetEntriesFast();
610 for (Int_t ientry=0; ientry<entries; ientry++) {
611 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
612 if (array) array->Delete();
613 delete fTrackHypothesys.RemoveAt(ientry);
616 fTrackHypothesys.Delete();
617 fBestHypothesys.Delete();
619 delete [] fCoefficients;
621 DeleteTrksMaterialLUT();
623 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
625 fTrackingPhase="Default";
629 //------------------------------------------------------------------------
630 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
631 //--------------------------------------------------------------------
632 // This functions propagates reconstructed ITS tracks back
633 // The clusters must be loaded !
634 //--------------------------------------------------------------------
635 fTrackingPhase="PropagateBack";
636 Int_t nentr=event->GetNumberOfTracks();
637 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
640 for (Int_t i=0; i<nentr; i++) {
641 AliESDtrack *esd=event->GetTrack(i);
643 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
644 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
648 t=new AliITStrackMI(*esd);
649 } catch (const Char_t *msg) {
650 //Warning("PropagateBack",msg);
654 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
656 ResetTrackToFollow(*t);
658 // propagate to vertex [SR, GSI 17.02.2003]
659 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
660 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
661 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
662 fTrackToFollow.StartTimeIntegral();
663 // from vertex to outside pipe
664 CorrectForPipeMaterial(&fTrackToFollow,"outward");
667 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
668 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
669 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
673 fTrackToFollow.SetLabel(t->GetLabel());
674 //fTrackToFollow.CookdEdx();
675 CookLabel(&fTrackToFollow,0.); //For comparison only
676 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
677 //UseClusters(&fTrackToFollow);
683 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
685 fTrackingPhase="Default";
689 //------------------------------------------------------------------------
690 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
691 //--------------------------------------------------------------------
692 // This functions refits ITS tracks using the
693 // "inward propagated" TPC tracks
694 // The clusters must be loaded !
695 //--------------------------------------------------------------------
696 fTrackingPhase="RefitInward";
697 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) RefitV02(event);
698 Int_t nentr=event->GetNumberOfTracks();
699 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
702 for (Int_t i=0; i<nentr; i++) {
703 AliESDtrack *esd=event->GetTrack(i);
705 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
706 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
707 if (esd->GetStatus()&AliESDtrack::kTPCout)
708 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
712 t=new AliITStrackMI(*esd);
713 } catch (const Char_t *msg) {
714 //Warning("RefitInward",msg);
718 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
719 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
724 ResetTrackToFollow(*t);
725 fTrackToFollow.ResetClusters();
727 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
728 fTrackToFollow.ResetCovariance(10.);
731 Bool_t pe=AliITSReconstructor::GetRecoParam()->GetComputePlaneEff();
732 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
733 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
734 AliDebug(2," refit OK");
735 fTrackToFollow.SetLabel(t->GetLabel());
736 // fTrackToFollow.CookdEdx();
737 CookdEdx(&fTrackToFollow);
739 CookLabel(&fTrackToFollow,0.0); //For comparison only
742 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
743 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
744 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
745 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
746 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
747 Float_t r[3]={0.,0.,0.};
749 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
756 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
758 fTrackingPhase="Default";
762 //------------------------------------------------------------------------
763 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
764 //--------------------------------------------------------------------
765 // Return pointer to a given cluster
766 //--------------------------------------------------------------------
767 Int_t l=(index & 0xf0000000) >> 28;
768 Int_t c=(index & 0x0fffffff) >> 00;
769 return fgLayers[l].GetCluster(c);
771 //------------------------------------------------------------------------
772 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
773 //--------------------------------------------------------------------
774 // Get track space point with index i
775 //--------------------------------------------------------------------
777 Int_t l=(index & 0xf0000000) >> 28;
778 Int_t c=(index & 0x0fffffff) >> 00;
779 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
780 Int_t idet = cl->GetDetectorIndex();
784 cl->GetGlobalXYZ(xyz);
785 cl->GetGlobalCov(cov);
787 p.SetCharge(cl->GetQ());
788 p.SetDriftTime(cl->GetDriftTime());
789 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
792 iLayer = AliGeomManager::kSPD1;
795 iLayer = AliGeomManager::kSPD2;
798 iLayer = AliGeomManager::kSDD1;
801 iLayer = AliGeomManager::kSDD2;
804 iLayer = AliGeomManager::kSSD1;
807 iLayer = AliGeomManager::kSSD2;
810 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
813 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
814 p.SetVolumeID((UShort_t)volid);
817 //------------------------------------------------------------------------
818 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
819 AliTrackPoint& p, const AliESDtrack *t) {
820 //--------------------------------------------------------------------
821 // Get track space point with index i
822 // (assign error estimated during the tracking)
823 //--------------------------------------------------------------------
825 Int_t l=(index & 0xf0000000) >> 28;
826 Int_t c=(index & 0x0fffffff) >> 00;
827 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
828 Int_t idet = cl->GetDetectorIndex();
830 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
832 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
834 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
835 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
836 Double_t alpha = t->GetAlpha();
837 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
838 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
839 phi += alpha-det.GetPhi();
840 Float_t tgphi = TMath::Tan(phi);
842 Float_t tgl = t->GetTgl(); // tgl about const along track
843 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
845 Float_t errlocalx,errlocalz;
846 Bool_t addMisalErr=kFALSE;
847 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
851 cl->GetGlobalXYZ(xyz);
852 // cl->GetGlobalCov(cov);
853 Float_t pos[3] = {0.,0.,0.};
854 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
855 tmpcl.GetGlobalCov(cov);
858 p.SetCharge(cl->GetQ());
859 p.SetDriftTime(cl->GetDriftTime());
861 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
864 iLayer = AliGeomManager::kSPD1;
867 iLayer = AliGeomManager::kSPD2;
870 iLayer = AliGeomManager::kSDD1;
873 iLayer = AliGeomManager::kSDD2;
876 iLayer = AliGeomManager::kSSD1;
879 iLayer = AliGeomManager::kSSD2;
882 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
885 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
887 p.SetVolumeID((UShort_t)volid);
890 //------------------------------------------------------------------------
891 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
893 //--------------------------------------------------------------------
894 // Follow prolongation tree
895 //--------------------------------------------------------------------
897 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
898 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
901 AliESDtrack * esd = otrack->GetESDtrack();
902 if (esd->GetV0Index(0)>0) {
903 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
904 // mapping of ESD track is different as ITS track in Containers
905 // Need something more stable
906 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
907 for (Int_t i=0;i<3;i++){
908 Int_t index = esd->GetV0Index(i);
910 AliESDv0 * vertex = fEsd->GetV0(index);
911 if (vertex->GetStatus()<0) continue; // rejected V0
913 if (esd->GetSign()>0) {
914 vertex->SetIndex(0,esdindex);
916 vertex->SetIndex(1,esdindex);
920 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
922 bestarray = new TObjArray(5);
923 fBestHypothesys.AddAt(bestarray,esdindex);
927 //setup tree of the prolongations
929 static AliITStrackMI tracks[7][100];
930 AliITStrackMI *currenttrack;
931 static AliITStrackMI currenttrack1;
932 static AliITStrackMI currenttrack2;
933 static AliITStrackMI backuptrack;
935 Int_t nindexes[7][100];
936 Float_t normalizedchi2[100];
937 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
938 otrack->SetNSkipped(0);
939 new (&(tracks[6][0])) AliITStrackMI(*otrack);
941 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
942 Int_t modstatus = 1; // found
946 // follow prolongations
947 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
948 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
951 AliITSlayer &layer=fgLayers[ilayer];
952 Double_t r = layer.GetR();
958 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
960 if (ntracks[ilayer]>=100) break;
961 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
962 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
963 if (ntracks[ilayer]>15+ilayer){
964 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
965 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
968 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
970 // material between SSD and SDD, SDD and SPD
972 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
974 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
978 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
979 Int_t idet=layer.FindDetectorIndex(phi,z);
981 Double_t trackGlobXYZ1[3];
982 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
984 // Get the budget to the primary vertex for the current track being prolonged
985 Double_t budgetToPrimVertex = GetEffectiveThickness();
987 // check if we allow a prolongation without point
988 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
990 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
991 // propagate to the layer radius
992 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
993 if(!vtrack->Propagate(xToGo)) continue;
994 // apply correction for material of the current layer
995 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
996 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
997 vtrack->SetClIndex(ilayer,0);
998 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
999 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1000 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1002 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1007 // track outside layer acceptance in z
1008 if (idet<0) continue;
1010 //propagate to the intersection with the detector plane
1011 const AliITSdetector &det=layer.GetDetector(idet);
1012 new(¤ttrack2) AliITStrackMI(currenttrack1);
1013 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1014 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1015 currenttrack1.SetDetectorIndex(idet);
1016 currenttrack2.SetDetectorIndex(idet);
1017 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1020 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1022 // road in global (rphi,z) [i.e. in tracking ref. system]
1023 Double_t zmin,zmax,ymin,ymax;
1024 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1026 // select clusters in road
1027 layer.SelectClusters(zmin,zmax,ymin,ymax);
1028 //********************
1030 // Define criteria for track-cluster association
1031 Double_t msz = currenttrack1.GetSigmaZ2() +
1032 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1033 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1034 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1035 Double_t msy = currenttrack1.GetSigmaY2() +
1036 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1037 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1038 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1040 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1041 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1043 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1044 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1046 msz = 1./msz; // 1/RoadZ^2
1047 msy = 1./msy; // 1/RoadY^2
1051 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1053 const AliITSRecPoint *cl=0;
1055 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1056 Bool_t deadzoneSPD=kFALSE;
1057 currenttrack = ¤ttrack1;
1059 // check if the road contains a dead zone
1060 Bool_t noClusters = kFALSE;
1061 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1062 if (noClusters) AliDebug(2,"no clusters in road");
1063 Double_t dz=0.5*(zmax-zmin);
1064 Double_t dy=0.5*(ymax-ymin);
1065 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1066 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1067 // create a prolongation without clusters (check also if there are no clusters in the road)
1070 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1071 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1072 updatetrack->SetClIndex(ilayer,0);
1074 modstatus = 5; // no cls in road
1075 } else if (dead==1) {
1076 modstatus = 7; // holes in z in SPD
1077 } else if (dead==2 || dead==3) {
1078 modstatus = 2; // dead from OCDB
1080 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1081 // apply correction for material of the current layer
1082 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1083 if (constrain) { // apply vertex constrain
1084 updatetrack->SetConstrain(constrain);
1085 Bool_t isPrim = kTRUE;
1086 if (ilayer<4) { // check that it's close to the vertex
1087 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1088 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1089 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1090 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1091 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1093 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1096 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1097 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1098 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1106 // loop over clusters in the road
1107 while ((cl=layer.GetNextCluster(clidx))!=0) {
1108 if (ntracks[ilayer]>95) break; //space for skipped clusters
1109 Bool_t changedet =kFALSE;
1110 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1111 Int_t idetc=cl->GetDetectorIndex();
1113 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1114 // take into account misalignment (bring track to real detector plane)
1115 Double_t xTrOrig = currenttrack->GetX();
1116 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1117 // a first cut on track-cluster distance
1118 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1119 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1120 { // cluster not associated to track
1121 AliDebug(2,"not associated");
1124 // bring track back to ideal detector plane
1125 if (!currenttrack->Propagate(xTrOrig)) continue;
1126 } else { // have to move track to cluster's detector
1127 const AliITSdetector &detc=layer.GetDetector(idetc);
1128 // a first cut on track-cluster distance
1130 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1131 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1132 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1133 continue; // cluster not associated to track
1135 new (&backuptrack) AliITStrackMI(currenttrack2);
1137 currenttrack =¤ttrack2;
1138 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1139 new (currenttrack) AliITStrackMI(backuptrack);
1143 currenttrack->SetDetectorIndex(idetc);
1144 // Get again the budget to the primary vertex
1145 // for the current track being prolonged, if had to change detector
1146 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1149 // calculate track-clusters chi2
1150 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1152 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1153 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1154 if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1155 if (ntracks[ilayer]>=100) continue;
1156 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1157 updatetrack->SetClIndex(ilayer,0);
1158 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1160 if (cl->GetQ()!=0) { // real cluster
1161 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1162 AliDebug(2,"update failed");
1165 updatetrack->SetSampledEdx(cl->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
1166 modstatus = 1; // found
1167 } else { // virtual cluster in dead zone
1168 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1169 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1170 modstatus = 7; // holes in z in SPD
1174 Float_t xlocnewdet,zlocnewdet;
1175 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1176 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1179 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1181 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1183 // apply correction for material of the current layer
1184 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1186 if (constrain) { // apply vertex constrain
1187 updatetrack->SetConstrain(constrain);
1188 Bool_t isPrim = kTRUE;
1189 if (ilayer<4) { // check that it's close to the vertex
1190 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1191 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1192 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1193 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1194 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1196 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1197 } //apply vertex constrain
1199 } // create new hypothesis
1201 AliDebug(2,"chi2 too large");
1204 } // loop over possible prolongations
1206 // allow one prolongation without clusters
1207 if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1208 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1209 // apply correction for material of the current layer
1210 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1211 vtrack->SetClIndex(ilayer,0);
1212 modstatus = 3; // skipped
1213 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1214 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1215 vtrack->IncrementNSkipped();
1219 // allow one prolongation without clusters for tracks with |tgl|>1.1
1220 if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) { //big theta - for low flux
1221 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1222 // apply correction for material of the current layer
1223 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1224 vtrack->SetClIndex(ilayer,0);
1225 modstatus = 3; // skipped
1226 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1227 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1228 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1233 } // loop over tracks in layer ilayer+1
1235 //loop over track candidates for the current layer
1241 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1242 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1243 if (normalizedchi2[itrack] <
1244 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1248 if (constrain) { // constrain
1249 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1251 } else { // !constrain
1252 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1257 // sort tracks by increasing normalized chi2
1258 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1259 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1260 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1261 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1262 } // end loop over layers
1266 // Now select tracks to be kept
1268 Int_t max = constrain ? 20 : 5;
1270 // tracks that reach layer 0 (SPD inner)
1271 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1272 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1273 if (track.GetNumberOfClusters()<2) continue;
1274 if (!constrain && track.GetNormChi2(0) >
1275 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1278 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1281 // tracks that reach layer 1 (SPD outer)
1282 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1283 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1284 if (track.GetNumberOfClusters()<4) continue;
1285 if (!constrain && track.GetNormChi2(1) >
1286 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1287 if (constrain) track.IncrementNSkipped();
1289 track.SetD(0,track.GetD(GetX(),GetY()));
1290 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1291 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1292 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1295 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1298 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1300 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1301 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1302 if (track.GetNumberOfClusters()<3) continue;
1303 if (!constrain && track.GetNormChi2(2) >
1304 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1305 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1307 track.SetD(0,track.GetD(GetX(),GetY()));
1308 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1309 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1310 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1313 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1319 // register best track of each layer - important for V0 finder
1321 for (Int_t ilayer=0;ilayer<5;ilayer++){
1322 if (ntracks[ilayer]==0) continue;
1323 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1324 if (track.GetNumberOfClusters()<1) continue;
1325 CookLabel(&track,0);
1326 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1330 // update TPC V0 information
1332 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1333 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1334 for (Int_t i=0;i<3;i++){
1335 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1336 if (index==0) break;
1337 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1338 if (vertex->GetStatus()<0) continue; // rejected V0
1340 if (otrack->GetSign()>0) {
1341 vertex->SetIndex(0,esdindex);
1344 vertex->SetIndex(1,esdindex);
1346 //find nearest layer with track info
1347 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1348 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1349 Int_t nearest = nearestold;
1350 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1351 if (ntracks[nearest]==0){
1356 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1357 if (nearestold<5&&nearest<5){
1358 Bool_t accept = track.GetNormChi2(nearest)<10;
1360 if (track.GetSign()>0) {
1361 vertex->SetParamP(track);
1362 vertex->Update(fprimvertex);
1363 //vertex->SetIndex(0,track.fESDtrack->GetID());
1364 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1366 vertex->SetParamN(track);
1367 vertex->Update(fprimvertex);
1368 //vertex->SetIndex(1,track.fESDtrack->GetID());
1369 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1371 vertex->SetStatus(vertex->GetStatus()+1);
1373 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1380 //------------------------------------------------------------------------
1381 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1383 //--------------------------------------------------------------------
1386 return fgLayers[layer];
1388 //------------------------------------------------------------------------
1389 AliITStrackerMI::AliITSlayer::AliITSlayer():
1414 //--------------------------------------------------------------------
1415 //default AliITSlayer constructor
1416 //--------------------------------------------------------------------
1417 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1418 fClusterWeight[i]=0;
1419 fClusterTracks[0][i]=-1;
1420 fClusterTracks[1][i]=-1;
1421 fClusterTracks[2][i]=-1;
1422 fClusterTracks[3][i]=-1;
1425 //------------------------------------------------------------------------
1426 AliITStrackerMI::AliITSlayer::
1427 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1452 //--------------------------------------------------------------------
1453 //main AliITSlayer constructor
1454 //--------------------------------------------------------------------
1455 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1456 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1458 //------------------------------------------------------------------------
1459 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1461 fPhiOffset(layer.fPhiOffset),
1462 fNladders(layer.fNladders),
1463 fZOffset(layer.fZOffset),
1464 fNdetectors(layer.fNdetectors),
1465 fDetectors(layer.fDetectors),
1470 fClustersCs(layer.fClustersCs),
1471 fClusterIndexCs(layer.fClusterIndexCs),
1475 fCurrentSlice(layer.fCurrentSlice),
1482 fAccepted(layer.fAccepted),
1486 //------------------------------------------------------------------------
1487 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1488 //--------------------------------------------------------------------
1489 // AliITSlayer destructor
1490 //--------------------------------------------------------------------
1491 delete [] fDetectors;
1492 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1493 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1494 fClusterWeight[i]=0;
1495 fClusterTracks[0][i]=-1;
1496 fClusterTracks[1][i]=-1;
1497 fClusterTracks[2][i]=-1;
1498 fClusterTracks[3][i]=-1;
1501 //------------------------------------------------------------------------
1502 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1503 //--------------------------------------------------------------------
1504 // This function removes loaded clusters
1505 //--------------------------------------------------------------------
1506 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1507 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1508 fClusterWeight[i]=0;
1509 fClusterTracks[0][i]=-1;
1510 fClusterTracks[1][i]=-1;
1511 fClusterTracks[2][i]=-1;
1512 fClusterTracks[3][i]=-1;
1518 //------------------------------------------------------------------------
1519 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1520 //--------------------------------------------------------------------
1521 // This function reset weights of the clusters
1522 //--------------------------------------------------------------------
1523 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1524 fClusterWeight[i]=0;
1525 fClusterTracks[0][i]=-1;
1526 fClusterTracks[1][i]=-1;
1527 fClusterTracks[2][i]=-1;
1528 fClusterTracks[3][i]=-1;
1530 for (Int_t i=0; i<fN;i++) {
1531 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1532 if (cl&&cl->IsUsed()) cl->Use();
1536 //------------------------------------------------------------------------
1537 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1538 //--------------------------------------------------------------------
1539 // This function calculates the road defined by the cluster density
1540 //--------------------------------------------------------------------
1542 for (Int_t i=0; i<fN; i++) {
1543 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1545 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1547 //------------------------------------------------------------------------
1548 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1549 //--------------------------------------------------------------------
1550 //This function adds a cluster to this layer
1551 //--------------------------------------------------------------------
1552 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1553 ::Error("InsertCluster","Too many clusters !\n");
1559 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1560 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1561 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1562 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1563 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1567 //------------------------------------------------------------------------
1568 void AliITStrackerMI::AliITSlayer::SortClusters()
1573 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1574 Float_t *z = new Float_t[fN];
1575 Int_t * index = new Int_t[fN];
1577 for (Int_t i=0;i<fN;i++){
1578 z[i] = fClusters[i]->GetZ();
1580 TMath::Sort(fN,z,index,kFALSE);
1581 for (Int_t i=0;i<fN;i++){
1582 clusters[i] = fClusters[index[i]];
1585 for (Int_t i=0;i<fN;i++){
1586 fClusters[i] = clusters[i];
1587 fZ[i] = fClusters[i]->GetZ();
1588 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1589 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1590 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1600 for (Int_t i=0;i<fN;i++){
1601 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1602 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1603 fClusterIndex[i] = i;
1607 fDy5 = (fYB[1]-fYB[0])/5.;
1608 fDy10 = (fYB[1]-fYB[0])/10.;
1609 fDy20 = (fYB[1]-fYB[0])/20.;
1610 for (Int_t i=0;i<6;i++) fN5[i] =0;
1611 for (Int_t i=0;i<11;i++) fN10[i]=0;
1612 for (Int_t i=0;i<21;i++) fN20[i]=0;
1614 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;}
1615 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;}
1616 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;}
1619 for (Int_t i=0;i<fN;i++)
1620 for (Int_t irot=-1;irot<=1;irot++){
1621 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1623 for (Int_t slice=0; slice<6;slice++){
1624 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1625 fClusters5[slice][fN5[slice]] = fClusters[i];
1626 fY5[slice][fN5[slice]] = curY;
1627 fZ5[slice][fN5[slice]] = fZ[i];
1628 fClusterIndex5[slice][fN5[slice]]=i;
1633 for (Int_t slice=0; slice<11;slice++){
1634 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1635 fClusters10[slice][fN10[slice]] = fClusters[i];
1636 fY10[slice][fN10[slice]] = curY;
1637 fZ10[slice][fN10[slice]] = fZ[i];
1638 fClusterIndex10[slice][fN10[slice]]=i;
1643 for (Int_t slice=0; slice<21;slice++){
1644 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1645 fClusters20[slice][fN20[slice]] = fClusters[i];
1646 fY20[slice][fN20[slice]] = curY;
1647 fZ20[slice][fN20[slice]] = fZ[i];
1648 fClusterIndex20[slice][fN20[slice]]=i;
1655 // consistency check
1657 for (Int_t i=0;i<fN-1;i++){
1663 for (Int_t slice=0;slice<21;slice++)
1664 for (Int_t i=0;i<fN20[slice]-1;i++){
1665 if (fZ20[slice][i]>fZ20[slice][i+1]){
1672 //------------------------------------------------------------------------
1673 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1674 //--------------------------------------------------------------------
1675 // This function returns the index of the nearest cluster
1676 //--------------------------------------------------------------------
1679 if (fCurrentSlice<0) {
1688 if (ncl==0) return 0;
1689 Int_t b=0, e=ncl-1, m=(b+e)/2;
1690 for (; b<e; m=(b+e)/2) {
1691 // if (z > fClusters[m]->GetZ()) b=m+1;
1692 if (z > zcl[m]) b=m+1;
1697 //------------------------------------------------------------------------
1698 Bool_t AliITStrackerMI::ComputeRoad(AliITStrackMI* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1699 //--------------------------------------------------------------------
1700 // This function computes the rectangular road for this track
1701 //--------------------------------------------------------------------
1704 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1705 // take into account the misalignment: propagate track to misaligned detector plane
1706 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1708 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1709 TMath::Sqrt(track->GetSigmaZ2() +
1710 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1711 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1712 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1713 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1714 TMath::Sqrt(track->GetSigmaY2() +
1715 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1716 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1717 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1719 // track at boundary between detectors, enlarge road
1720 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1721 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1722 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1723 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1724 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1725 Float_t tgl = TMath::Abs(track->GetTgl());
1726 if (tgl > 1.) tgl=1.;
1727 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1728 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1729 Float_t snp = TMath::Abs(track->GetSnp());
1730 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1731 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1734 // add to the road a term (up to 2-3 mm) to deal with misalignments
1735 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1736 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1738 Double_t r = fgLayers[ilayer].GetR();
1739 zmin = track->GetZ() - dz;
1740 zmax = track->GetZ() + dz;
1741 ymin = track->GetY() + r*det.GetPhi() - dy;
1742 ymax = track->GetY() + r*det.GetPhi() + dy;
1744 // bring track back to idead detector plane
1745 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1749 //------------------------------------------------------------------------
1750 void AliITStrackerMI::AliITSlayer::
1751 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1752 //--------------------------------------------------------------------
1753 // This function sets the "window"
1754 //--------------------------------------------------------------------
1756 Double_t circle=2*TMath::Pi()*fR;
1757 fYmin = ymin; fYmax =ymax;
1758 Float_t ymiddle = (fYmin+fYmax)*0.5;
1759 if (ymiddle<fYB[0]) {
1760 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1761 } else if (ymiddle>fYB[1]) {
1762 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1768 fClustersCs = fClusters;
1769 fClusterIndexCs = fClusterIndex;
1775 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1776 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1777 if (slice<0) slice=0;
1778 if (slice>20) slice=20;
1779 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1781 fCurrentSlice=slice;
1782 fClustersCs = fClusters20[fCurrentSlice];
1783 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1784 fYcs = fY20[fCurrentSlice];
1785 fZcs = fZ20[fCurrentSlice];
1786 fNcs = fN20[fCurrentSlice];
1791 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1792 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1793 if (slice<0) slice=0;
1794 if (slice>10) slice=10;
1795 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1797 fCurrentSlice=slice;
1798 fClustersCs = fClusters10[fCurrentSlice];
1799 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1800 fYcs = fY10[fCurrentSlice];
1801 fZcs = fZ10[fCurrentSlice];
1802 fNcs = fN10[fCurrentSlice];
1807 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1808 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1809 if (slice<0) slice=0;
1810 if (slice>5) slice=5;
1811 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1813 fCurrentSlice=slice;
1814 fClustersCs = fClusters5[fCurrentSlice];
1815 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1816 fYcs = fY5[fCurrentSlice];
1817 fZcs = fZ5[fCurrentSlice];
1818 fNcs = fN5[fCurrentSlice];
1822 fI=FindClusterIndex(zmin); fZmax=zmax;
1823 fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1829 //------------------------------------------------------------------------
1830 Int_t AliITStrackerMI::AliITSlayer::
1831 FindDetectorIndex(Double_t phi, Double_t z) const {
1832 //--------------------------------------------------------------------
1833 //This function finds the detector crossed by the track
1834 //--------------------------------------------------------------------
1836 if (fZOffset<0) // old geometry
1837 dphi = -(phi-fPhiOffset);
1838 else // new geometry
1839 dphi = phi-fPhiOffset;
1842 if (dphi < 0) dphi += 2*TMath::Pi();
1843 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1844 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1845 if (np>=fNladders) np-=fNladders;
1846 if (np<0) np+=fNladders;
1849 Double_t dz=fZOffset-z;
1850 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1851 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1852 if (nz>=fNdetectors) return -1;
1853 if (nz<0) return -1;
1855 // ad hoc correction for 3rd ladder of SDD inner layer,
1856 // which is reversed (rotated by pi around local y)
1857 // this correction is OK only from AliITSv11Hybrid onwards
1858 if (GetR()>12. && GetR()<20.) { // SDD inner
1859 if(np==2) { // 3rd ladder
1860 nz = (fNdetectors-1) - nz;
1863 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1866 return np*fNdetectors + nz;
1868 //------------------------------------------------------------------------
1869 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1871 //--------------------------------------------------------------------
1872 // This function returns clusters within the "window"
1873 //--------------------------------------------------------------------
1875 if (fCurrentSlice<0) {
1876 Double_t rpi2 = 2.*fR*TMath::Pi();
1877 for (Int_t i=fI; i<fImax; i++) {
1879 if (fYmax<y) y -= rpi2;
1880 if (fYmin>y) y += rpi2;
1881 if (y<fYmin) continue;
1882 if (y>fYmax) continue;
1883 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1886 return fClusters[i];
1889 for (Int_t i=fI; i<fImax; i++) {
1890 if (fYcs[i]<fYmin) continue;
1891 if (fYcs[i]>fYmax) continue;
1892 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1893 ci=fClusterIndexCs[i];
1895 return fClustersCs[i];
1900 //------------------------------------------------------------------------
1901 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1903 //--------------------------------------------------------------------
1904 // This function returns the layer thickness at this point (units X0)
1905 //--------------------------------------------------------------------
1907 x0=AliITSRecoParam::GetX0Air();
1908 if (43<fR&&fR<45) { //SSD2
1911 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1912 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1913 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1914 for (Int_t i=0; i<12; i++) {
1915 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1916 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1920 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1921 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1925 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1926 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1929 if (37<fR&&fR<41) { //SSD1
1932 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1933 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1934 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1935 for (Int_t i=0; i<11; i++) {
1936 if (TMath::Abs(z-3.9*i)<0.15) {
1937 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1941 if (TMath::Abs(z+3.9*i)<0.15) {
1942 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1946 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1947 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1950 if (13<fR&&fR<26) { //SDD
1953 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1955 if (TMath::Abs(y-1.80)<0.55) {
1957 for (Int_t j=0; j<20; j++) {
1958 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1959 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1962 if (TMath::Abs(y+1.80)<0.55) {
1964 for (Int_t j=0; j<20; j++) {
1965 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1966 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1970 for (Int_t i=0; i<4; i++) {
1971 if (TMath::Abs(z-7.3*i)<0.60) {
1973 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1976 if (TMath::Abs(z+7.3*i)<0.60) {
1978 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1983 if (6<fR&&fR<8) { //SPD2
1984 Double_t dd=0.0063; x0=21.5;
1986 if (TMath::Abs(y-3.08)>0.5) d+=dd;
1987 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
1989 if (3<fR&&fR<5) { //SPD1
1990 Double_t dd=0.0063; x0=21.5;
1992 if (TMath::Abs(y+0.21)>0.6) d+=dd;
1993 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
1998 //------------------------------------------------------------------------
1999 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2001 fRmisal(det.fRmisal),
2003 fSinPhi(det.fSinPhi),
2004 fCosPhi(det.fCosPhi),
2010 fNChips(det.fNChips),
2011 fChipIsBad(det.fChipIsBad)
2015 //------------------------------------------------------------------------
2016 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2017 AliITSDetTypeRec *detTypeRec)
2019 //--------------------------------------------------------------------
2020 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2021 //--------------------------------------------------------------------
2023 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2024 // while in the tracker they start from 0 for each layer
2025 for(Int_t il=0; il<ilayer; il++)
2026 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2029 if (ilayer==0 || ilayer==1) { // ---------- SPD
2031 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2033 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2036 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2040 // Get calibration from AliITSDetTypeRec
2041 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2042 calib->SetModuleIndex(idet);
2043 AliITSCalibration *calibSPDdead = 0;
2044 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2045 if (calib->IsBad() ||
2046 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2049 // printf("lay %d bad %d\n",ilayer,idet);
2052 // Get segmentation from AliITSDetTypeRec
2053 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2055 // Read info about bad chips
2056 fNChips = segm->GetMaximumChipIndex()+1;
2057 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2058 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2059 fChipIsBad = new Bool_t[fNChips];
2060 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2061 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2062 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2063 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2068 //------------------------------------------------------------------------
2069 Double_t AliITStrackerMI::GetEffectiveThickness()
2071 //--------------------------------------------------------------------
2072 // Returns the thickness between the current layer and the vertex (units X0)
2073 //--------------------------------------------------------------------
2076 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2077 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2078 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2082 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2083 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2087 Double_t xn=fgLayers[fI].GetR();
2088 for (Int_t i=0; i<fI; i++) {
2089 Double_t xi=fgLayers[i].GetR();
2090 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2096 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2097 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2100 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2101 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2105 //------------------------------------------------------------------------
2106 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2107 //-------------------------------------------------------------------
2108 // This function returns number of clusters within the "window"
2109 //--------------------------------------------------------------------
2111 for (Int_t i=fI; i<fN; i++) {
2112 const AliITSRecPoint *c=fClusters[i];
2113 if (c->GetZ() > fZmax) break;
2114 if (c->IsUsed()) continue;
2115 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2116 Double_t y=fR*det.GetPhi() + c->GetY();
2118 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2119 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2121 if (y<fYmin) continue;
2122 if (y>fYmax) continue;
2127 //------------------------------------------------------------------------
2128 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2129 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2131 //--------------------------------------------------------------------
2132 // This function refits the track "track" at the position "x" using
2133 // the clusters from "clusters"
2134 // If "extra"==kTRUE,
2135 // the clusters from overlapped modules get attached to "track"
2136 // If "planeff"==kTRUE,
2137 // special approach for plane efficiency evaluation is applyed
2138 //--------------------------------------------------------------------
2140 Int_t index[AliITSgeomTGeo::kNLayers];
2142 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2143 Int_t nc=clusters->GetNumberOfClusters();
2144 for (k=0; k<nc; k++) {
2145 Int_t idx=clusters->GetClusterIndex(k);
2146 Int_t ilayer=(idx&0xf0000000)>>28;
2150 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2152 //------------------------------------------------------------------------
2153 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2154 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2156 //--------------------------------------------------------------------
2157 // This function refits the track "track" at the position "x" using
2158 // the clusters from array
2159 // If "extra"==kTRUE,
2160 // the clusters from overlapped modules get attached to "track"
2161 // If "planeff"==kTRUE,
2162 // special approach for plane efficiency evaluation is applyed
2163 //--------------------------------------------------------------------
2164 Int_t index[AliITSgeomTGeo::kNLayers];
2166 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2168 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2169 index[k]=clusters[k];
2172 // special for cosmics: check which the innermost layer crossed
2174 Int_t innermostlayer=5;
2175 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2176 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2177 if(drphi < fgLayers[innermostlayer].GetR()) break;
2179 //printf(" drphi %f innermost %d\n",drphi,innermostlayer);
2181 Int_t modstatus=1; // found
2183 Int_t from, to, step;
2184 if (xx > track->GetX()) {
2185 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2188 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2191 TString dir = (step>0 ? "outward" : "inward");
2193 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2194 AliITSlayer &layer=fgLayers[ilayer];
2195 Double_t r=layer.GetR();
2196 if (step<0 && xx>r) break;
2198 // material between SSD and SDD, SDD and SPD
2199 Double_t hI=ilayer-0.5*step;
2200 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2201 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2202 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2203 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2205 // remember old position [SR, GSI 18.02.2003]
2206 Double_t oldX=0., oldY=0., oldZ=0.;
2207 if (track->IsStartedTimeIntegral() && step==1) {
2208 if (!track->GetGlobalXYZat(track->GetX(),oldX,oldY,oldZ)) return kFALSE;
2212 Double_t oldGlobXYZ[3];
2213 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2214 //TMath::Sqrt(track->GetSigmaY2());
2217 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2219 Int_t idet=layer.FindDetectorIndex(phi,z);
2221 // check if we allow a prolongation without point for large-eta tracks
2222 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2224 // propagate to the layer radius
2225 Double_t xToGo; if (!track->GetLocalXat(r,xToGo)) return kFALSE;
2226 if (!track->Propagate(xToGo)) return kFALSE;
2227 // apply correction for material of the current layer
2228 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2229 modstatus = 4; // out in z
2230 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2231 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2233 // track time update [SR, GSI 17.02.2003]
2234 if (track->IsStartedTimeIntegral() && step==1) {
2235 Double_t newX, newY, newZ;
2236 if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2237 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2238 (oldZ-newZ)*(oldZ-newZ);
2239 track->AddTimeStep(TMath::Sqrt(dL2));
2244 if (idet<0) return kFALSE;
2246 const AliITSdetector &det=layer.GetDetector(idet);
2247 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2249 track->SetDetectorIndex(idet);
2250 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2252 Double_t dz,zmin,zmax,dy,ymin,ymax;
2254 const AliITSRecPoint *clAcc=0;
2255 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2257 Int_t idx=index[ilayer];
2258 if (idx>=0) { // cluster in this layer
2259 modstatus = 6; // no refit
2260 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2262 if (idet != cl->GetDetectorIndex()) {
2263 idet=cl->GetDetectorIndex();
2264 const AliITSdetector &detc=layer.GetDetector(idet);
2265 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2266 track->SetDetectorIndex(idet);
2267 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2269 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2270 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2274 modstatus = 1; // found
2279 } else { // no cluster in this layer
2281 modstatus = 3; // skipped
2282 // Plane Eff determination:
2283 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2284 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2285 UseTrackForPlaneEff(track,ilayer);
2288 modstatus = 5; // no cls in road
2290 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2291 dz = 0.5*(zmax-zmin);
2292 dy = 0.5*(ymax-ymin);
2293 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2294 if (dead==1) modstatus = 7; // holes in z in SPD
2295 if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2300 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2301 track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2303 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2306 if (extra) { // search for extra clusters in overlapped modules
2307 AliITStrackV2 tmp(*track);
2308 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2309 layer.SelectClusters(zmin,zmax,ymin,ymax);
2311 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2313 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2314 Double_t tolerance=0.1;
2315 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2316 // only clusters in another module! (overlaps)
2317 idetExtra = clExtra->GetDetectorIndex();
2318 if (idet == idetExtra) continue;
2320 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2322 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2323 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2324 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2325 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2327 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2328 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2331 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2332 track->SetExtraModule(ilayer,idetExtra);
2334 } // end search for extra clusters in overlapped modules
2336 // Correct for material of the current layer
2337 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2339 // track time update [SR, GSI 17.02.2003]
2340 if (track->IsStartedTimeIntegral() && step==1) {
2341 Double_t newX, newY, newZ;
2342 if (!track->GetGlobalXYZat(track->GetX(),newX,newY,newZ)) return kFALSE;
2343 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2344 (oldZ-newZ)*(oldZ-newZ);
2345 track->AddTimeStep(TMath::Sqrt(dL2));
2349 } // end loop on layers
2351 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2355 //------------------------------------------------------------------------
2356 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2359 // calculate normalized chi2
2360 // return NormalizedChi2(track,0);
2363 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2364 // track->fdEdxMismatch=0;
2365 Float_t dedxmismatch =0;
2366 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2368 for (Int_t i = 0;i<6;i++){
2369 if (track->GetClIndex(i)>0){
2370 Float_t cerry, cerrz;
2371 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2373 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2376 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2377 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2378 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2380 cchi2+=(0.5-ratio)*10.;
2381 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2382 dedxmismatch+=(0.5-ratio)*10.;
2386 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2387 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2388 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2389 if (i<2) chi2+=2*cl->GetDeltaProbability();
2395 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2396 track->SetdEdxMismatch(dedxmismatch);
2400 for (Int_t i = 0;i<4;i++){
2401 if (track->GetClIndex(i)>0){
2402 Float_t cerry, cerrz;
2403 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2404 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2407 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2408 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2412 for (Int_t i = 4;i<6;i++){
2413 if (track->GetClIndex(i)>0){
2414 Float_t cerry, cerrz;
2415 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2416 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2419 Float_t cerryb, cerrzb;
2420 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2421 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2424 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2425 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2430 if (track->GetESDtrack()->GetTPCsignal()>85){
2431 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2433 chi2+=(0.5-ratio)*5.;
2436 chi2+=(ratio-2.0)*3;
2440 Double_t match = TMath::Sqrt(track->GetChi22());
2441 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2442 if (!track->GetConstrain()) {
2443 if (track->GetNumberOfClusters()>2) {
2444 match/=track->GetNumberOfClusters()-2.;
2449 if (match<0) match=0;
2450 Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2451 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2452 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2453 1./(1.+track->GetNSkipped()));
2457 //------------------------------------------------------------------------
2458 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2461 // return matching chi2 between two tracks
2462 Double_t largeChi2=1000.;
2464 AliITStrackMI track3(*track2);
2465 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2467 vec(0,0)=track1->GetY() - track3.GetY();
2468 vec(1,0)=track1->GetZ() - track3.GetZ();
2469 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2470 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2471 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2474 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2475 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2476 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2477 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2478 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2480 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2481 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2482 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2483 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2485 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2486 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2487 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2489 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2490 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2492 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2495 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2496 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2499 //------------------------------------------------------------------------
2500 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr)
2503 // return probability that given point (characterized by z position and error)
2504 // is in SPD dead zone
2506 Double_t probability = 0.;
2507 Double_t absz = TMath::Abs(zpos);
2508 Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2509 0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2510 if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2511 Double_t zmin, zmax;
2512 if (zpos<-6.) { // dead zone at z = -7
2513 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2514 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2515 } else if (zpos>6.) { // dead zone at z = +7
2516 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2517 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2518 } else if (absz<2.) { // dead zone at z = 0
2519 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2520 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2525 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2527 probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2528 TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2531 //------------------------------------------------------------------------
2532 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2535 // calculate normalized chi2
2537 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2539 for (Int_t i = 0;i<6;i++){
2540 if (TMath::Abs(track->GetDy(i))>0){
2541 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2542 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2545 else{chi2[i]=10000;}
2548 TMath::Sort(6,chi2,index,kFALSE);
2549 Float_t max = float(ncl)*fac-1.;
2550 Float_t sumchi=0, sumweight=0;
2551 for (Int_t i=0;i<max+1;i++){
2552 Float_t weight = (i<max)?1.:(max+1.-i);
2553 sumchi+=weight*chi2[index[i]];
2556 Double_t normchi2 = sumchi/sumweight;
2559 //------------------------------------------------------------------------
2560 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2563 // calculate normalized chi2
2564 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2567 for (Int_t i=0;i<6;i++){
2568 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2569 Double_t sy1 = forwardtrack->GetSigmaY(i);
2570 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2571 Double_t sy2 = backtrack->GetSigmaY(i);
2572 Double_t sz2 = backtrack->GetSigmaZ(i);
2573 if (i<2){ sy2=1000.;sz2=1000;}
2575 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2576 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2578 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2579 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2581 res+= nz0*nz0+ny0*ny0;
2584 if (npoints>1) return
2585 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2586 //2*forwardtrack->fNUsed+
2587 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2588 1./(1.+forwardtrack->GetNSkipped()));
2591 //------------------------------------------------------------------------
2592 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2593 //--------------------------------------------------------------------
2594 // Return pointer to a given cluster
2595 //--------------------------------------------------------------------
2596 Int_t l=(index & 0xf0000000) >> 28;
2597 Int_t c=(index & 0x0fffffff) >> 00;
2598 return fgLayers[l].GetWeight(c);
2600 //------------------------------------------------------------------------
2601 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2603 //---------------------------------------------
2604 // register track to the list
2606 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2609 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2610 Int_t index = track->GetClusterIndex(icluster);
2611 Int_t l=(index & 0xf0000000) >> 28;
2612 Int_t c=(index & 0x0fffffff) >> 00;
2613 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2614 for (Int_t itrack=0;itrack<4;itrack++){
2615 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2616 fgLayers[l].SetClusterTracks(itrack,c,id);
2622 //------------------------------------------------------------------------
2623 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2625 //---------------------------------------------
2626 // unregister track from the list
2627 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2628 Int_t index = track->GetClusterIndex(icluster);
2629 Int_t l=(index & 0xf0000000) >> 28;
2630 Int_t c=(index & 0x0fffffff) >> 00;
2631 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2632 for (Int_t itrack=0;itrack<4;itrack++){
2633 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2634 fgLayers[l].SetClusterTracks(itrack,c,-1);
2639 //------------------------------------------------------------------------
2640 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2642 //-------------------------------------------------------------
2643 //get number of shared clusters
2644 //-------------------------------------------------------------
2646 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2647 // mean number of clusters
2648 Float_t *ny = GetNy(id), *nz = GetNz(id);
2651 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2652 Int_t index = track->GetClusterIndex(icluster);
2653 Int_t l=(index & 0xf0000000) >> 28;
2654 Int_t c=(index & 0x0fffffff) >> 00;
2655 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2657 printf("problem\n");
2659 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2663 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2664 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2665 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2667 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2670 deltan = (cl->GetNz()-nz[l]);
2672 if (deltan>2.0) continue; // extended - highly probable shared cluster
2673 weight = 2./TMath::Max(3.+deltan,2.);
2675 for (Int_t itrack=0;itrack<4;itrack++){
2676 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2678 clist[l] = (AliITSRecPoint*)GetCluster(index);
2684 track->SetNUsed(shared);
2687 //------------------------------------------------------------------------
2688 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2691 // find first shared track
2693 // mean number of clusters
2694 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2696 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2697 Int_t sharedtrack=100000;
2698 Int_t tracks[24],trackindex=0;
2699 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2701 for (Int_t icluster=0;icluster<6;icluster++){
2702 if (clusterlist[icluster]<0) continue;
2703 Int_t index = clusterlist[icluster];
2704 Int_t l=(index & 0xf0000000) >> 28;
2705 Int_t c=(index & 0x0fffffff) >> 00;
2707 printf("problem\n");
2709 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2710 //if (l>3) continue;
2711 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2714 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2715 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2716 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2718 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2721 deltan = (cl->GetNz()-nz[l]);
2723 if (deltan>2.0) continue; // extended - highly probable shared cluster
2725 for (Int_t itrack=3;itrack>=0;itrack--){
2726 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2727 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2728 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2733 if (trackindex==0) return -1;
2735 sharedtrack = tracks[0];
2737 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2740 Int_t tracks2[24], cluster[24];
2741 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2744 for (Int_t i=0;i<trackindex;i++){
2745 if (tracks[i]<0) continue;
2746 tracks2[index] = tracks[i];
2748 for (Int_t j=i+1;j<trackindex;j++){
2749 if (tracks[j]<0) continue;
2750 if (tracks[j]==tracks[i]){
2758 for (Int_t i=0;i<index;i++){
2759 if (cluster[index]>max) {
2760 sharedtrack=tracks2[index];
2767 if (sharedtrack>=100000) return -1;
2769 // make list of overlaps
2771 for (Int_t icluster=0;icluster<6;icluster++){
2772 if (clusterlist[icluster]<0) continue;
2773 Int_t index = clusterlist[icluster];
2774 Int_t l=(index & 0xf0000000) >> 28;
2775 Int_t c=(index & 0x0fffffff) >> 00;
2776 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2777 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2779 if (cl->GetNy()>2) continue;
2780 if (cl->GetNz()>2) continue;
2783 if (cl->GetNy()>3) continue;
2784 if (cl->GetNz()>3) continue;
2787 for (Int_t itrack=3;itrack>=0;itrack--){
2788 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2789 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2797 //------------------------------------------------------------------------
2798 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2800 // try to find track hypothesys without conflicts
2801 // with minimal chi2;
2802 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2803 Int_t entries1 = arr1->GetEntriesFast();
2804 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2805 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2806 Int_t entries2 = arr2->GetEntriesFast();
2807 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2809 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2810 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2811 if (track10->Pt()>0.5+track20->Pt()) return track10;
2813 for (Int_t itrack=0;itrack<entries1;itrack++){
2814 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2815 UnRegisterClusterTracks(track,trackID1);
2818 for (Int_t itrack=0;itrack<entries2;itrack++){
2819 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2820 UnRegisterClusterTracks(track,trackID2);
2824 Float_t maxconflicts=6;
2825 Double_t maxchi2 =1000.;
2827 // get weight of hypothesys - using best hypothesy
2830 Int_t list1[6],list2[6];
2831 AliITSRecPoint *clist1[6], *clist2[6] ;
2832 RegisterClusterTracks(track10,trackID1);
2833 RegisterClusterTracks(track20,trackID2);
2834 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2835 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2836 UnRegisterClusterTracks(track10,trackID1);
2837 UnRegisterClusterTracks(track20,trackID2);
2840 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2841 Float_t nerry[6],nerrz[6];
2842 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2843 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2844 for (Int_t i=0;i<6;i++){
2845 if ( (erry1[i]>0) && (erry2[i]>0)) {
2846 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2847 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2849 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2850 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2852 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2853 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2854 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2857 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2858 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2859 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2867 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2868 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2869 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2870 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2872 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2873 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2874 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2876 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2877 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2878 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2881 Double_t sumw = w1+w2;
2885 w1 = (d2+0.5)/(d1+d2+1.);
2886 w2 = (d1+0.5)/(d1+d2+1.);
2888 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2889 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2891 // get pair of "best" hypothesys
2893 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2894 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2896 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2897 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2898 //if (track1->fFakeRatio>0) continue;
2899 RegisterClusterTracks(track1,trackID1);
2900 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2901 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2903 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2904 //if (track2->fFakeRatio>0) continue;
2906 RegisterClusterTracks(track2,trackID2);
2907 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2908 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2909 UnRegisterClusterTracks(track2,trackID2);
2911 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2912 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2913 if (nskipped>0.5) continue;
2915 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2916 if (conflict1+1<cconflict1) continue;
2917 if (conflict2+1<cconflict2) continue;
2921 for (Int_t i=0;i<6;i++){
2923 Float_t c1 =0.; // conflict coeficients
2925 if (clist1[i]&&clist2[i]){
2928 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2931 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2933 c1 = 2./TMath::Max(3.+deltan,2.);
2934 c2 = 2./TMath::Max(3.+deltan,2.);
2940 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2943 deltan = (clist1[i]->GetNz()-nz1[i]);
2945 c1 = 2./TMath::Max(3.+deltan,2.);
2952 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2955 deltan = (clist2[i]->GetNz()-nz2[i]);
2957 c2 = 2./TMath::Max(3.+deltan,2.);
2963 if (TMath::Abs(track1->GetDy(i))>0.) {
2964 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2965 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2966 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2967 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2969 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2972 if (TMath::Abs(track2->GetDy(i))>0.) {
2973 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2974 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2975 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2976 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2979 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2981 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2982 if (chi21>0) sum+=w1;
2983 if (chi22>0) sum+=w2;
2986 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2987 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2988 Double_t normchi2 = 2*conflict+sumchi2/sum;
2989 if ( normchi2 <maxchi2 ){
2992 maxconflicts = conflict;
2996 UnRegisterClusterTracks(track1,trackID1);
2999 // if (maxconflicts<4 && maxchi2<th0){
3000 if (maxchi2<th0*2.){
3001 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3002 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3003 track1->SetChi2MIP(5,maxconflicts);
3004 track1->SetChi2MIP(6,maxchi2);
3005 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3006 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3007 track1->SetChi2MIP(8,index1);
3008 fBestTrackIndex[trackID1] =index1;
3009 UpdateESDtrack(track1, AliESDtrack::kITSin);
3011 else if (track10->GetChi2MIP(0)<th1){
3012 track10->SetChi2MIP(5,maxconflicts);
3013 track10->SetChi2MIP(6,maxchi2);
3014 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3015 UpdateESDtrack(track10,AliESDtrack::kITSin);
3018 for (Int_t itrack=0;itrack<entries1;itrack++){
3019 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3020 UnRegisterClusterTracks(track,trackID1);
3023 for (Int_t itrack=0;itrack<entries2;itrack++){
3024 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3025 UnRegisterClusterTracks(track,trackID2);
3028 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3029 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3030 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3031 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3032 RegisterClusterTracks(track10,trackID1);
3034 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3035 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3036 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3037 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3038 RegisterClusterTracks(track20,trackID2);
3043 //------------------------------------------------------------------------
3044 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3045 //--------------------------------------------------------------------
3046 // This function marks clusters assigned to the track
3047 //--------------------------------------------------------------------
3048 AliTracker::UseClusters(t,from);
3050 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3051 //if (c->GetQ()>2) c->Use();
3052 if (c->GetSigmaZ2()>0.1) c->Use();
3053 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3054 //if (c->GetQ()>2) c->Use();
3055 if (c->GetSigmaZ2()>0.1) c->Use();
3058 //------------------------------------------------------------------------
3059 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3061 //------------------------------------------------------------------
3062 // add track to the list of hypothesys
3063 //------------------------------------------------------------------
3065 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3066 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3068 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3070 array = new TObjArray(10);
3071 fTrackHypothesys.AddAt(array,esdindex);
3073 array->AddLast(track);
3075 //------------------------------------------------------------------------
3076 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3078 //-------------------------------------------------------------------
3079 // compress array of track hypothesys
3080 // keep only maxsize best hypothesys
3081 //-------------------------------------------------------------------
3082 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3083 if (! (fTrackHypothesys.At(esdindex)) ) return;
3084 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3085 Int_t entries = array->GetEntriesFast();
3087 //- find preliminary besttrack as a reference
3088 Float_t minchi2=10000;
3090 AliITStrackMI * besttrack=0;
3091 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3092 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3093 if (!track) continue;
3094 Float_t chi2 = NormalizedChi2(track,0);
3096 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3097 track->SetLabel(tpcLabel);
3099 track->SetFakeRatio(1.);
3100 CookLabel(track,0.); //For comparison only
3102 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3103 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3104 if (track->GetNumberOfClusters()<maxn) continue;
3105 maxn = track->GetNumberOfClusters();
3112 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3113 delete array->RemoveAt(itrack);
3117 if (!besttrack) return;
3120 //take errors of best track as a reference
3121 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3122 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3123 for (Int_t j=0;j<6;j++) {
3124 if (besttrack->GetClIndex(j)>0){
3125 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3126 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3127 ny[j] = besttrack->GetNy(j);
3128 nz[j] = besttrack->GetNz(j);
3132 // calculate normalized chi2
3134 Float_t * chi2 = new Float_t[entries];
3135 Int_t * index = new Int_t[entries];
3136 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3137 for (Int_t itrack=0;itrack<entries;itrack++){
3138 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3140 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3141 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3142 chi2[itrack] = track->GetChi2MIP(0);
3144 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3145 delete array->RemoveAt(itrack);
3151 TMath::Sort(entries,chi2,index,kFALSE);
3152 besttrack = (AliITStrackMI*)array->At(index[0]);
3153 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3154 for (Int_t j=0;j<6;j++){
3155 if (besttrack->GetClIndex(j)>0){
3156 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3157 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3158 ny[j] = besttrack->GetNy(j);
3159 nz[j] = besttrack->GetNz(j);
3164 // calculate one more time with updated normalized errors
3165 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3166 for (Int_t itrack=0;itrack<entries;itrack++){
3167 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3169 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3170 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3171 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3174 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3175 delete array->RemoveAt(itrack);
3180 entries = array->GetEntriesFast();
3184 TObjArray * newarray = new TObjArray();
3185 TMath::Sort(entries,chi2,index,kFALSE);
3186 besttrack = (AliITStrackMI*)array->At(index[0]);
3189 for (Int_t j=0;j<6;j++){
3190 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3191 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3192 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3193 ny[j] = besttrack->GetNy(j);
3194 nz[j] = besttrack->GetNz(j);
3197 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3198 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3199 Float_t minn = besttrack->GetNumberOfClusters()-3;
3201 for (Int_t i=0;i<entries;i++){
3202 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3203 if (!track) continue;
3204 if (accepted>maxcut) break;
3205 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3206 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3207 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3208 delete array->RemoveAt(index[i]);
3212 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3213 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3214 if (!shortbest) accepted++;
3216 newarray->AddLast(array->RemoveAt(index[i]));
3217 for (Int_t j=0;j<6;j++){
3219 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3220 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3221 ny[j] = track->GetNy(j);
3222 nz[j] = track->GetNz(j);
3227 delete array->RemoveAt(index[i]);
3231 delete fTrackHypothesys.RemoveAt(esdindex);
3232 fTrackHypothesys.AddAt(newarray,esdindex);
3236 delete fTrackHypothesys.RemoveAt(esdindex);
3242 //------------------------------------------------------------------------
3243 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3245 //-------------------------------------------------------------
3246 // try to find best hypothesy
3247 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3248 //-------------------------------------------------------------
3249 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3250 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3251 if (!array) return 0;
3252 Int_t entries = array->GetEntriesFast();
3253 if (!entries) return 0;
3254 Float_t minchi2 = 100000;
3255 AliITStrackMI * besttrack=0;
3257 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3258 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3259 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3260 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3262 for (Int_t i=0;i<entries;i++){
3263 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3264 if (!track) continue;
3265 Float_t sigmarfi,sigmaz;
3266 GetDCASigma(track,sigmarfi,sigmaz);
3267 track->SetDnorm(0,sigmarfi);
3268 track->SetDnorm(1,sigmaz);
3270 track->SetChi2MIP(1,1000000);
3271 track->SetChi2MIP(2,1000000);
3272 track->SetChi2MIP(3,1000000);
3275 backtrack = new(backtrack) AliITStrackMI(*track);
3276 if (track->GetConstrain()) {
3277 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3278 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3279 backtrack->ResetCovariance(10.);
3281 backtrack->ResetCovariance(10.);
3283 backtrack->ResetClusters();
3285 Double_t x = original->GetX();
3286 if (!RefitAt(x,backtrack,track)) continue;
3288 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3289 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3290 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3291 track->SetChi22(GetMatchingChi2(backtrack,original));
3293 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3294 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3295 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3298 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3300 //forward track - without constraint
3301 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3302 forwardtrack->ResetClusters();
3304 RefitAt(x,forwardtrack,track);
3305 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3306 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3307 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3309 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3310 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3311 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3312 forwardtrack->SetD(0,track->GetD(0));
3313 forwardtrack->SetD(1,track->GetD(1));
3316 AliITSRecPoint* clist[6];
3317 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3318 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3321 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3322 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3323 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3324 track->SetChi2MIP(3,1000);
3327 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3329 for (Int_t ichi=0;ichi<5;ichi++){
3330 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3332 if (chi2 < minchi2){
3333 //besttrack = new AliITStrackMI(*forwardtrack);
3335 besttrack->SetLabel(track->GetLabel());
3336 besttrack->SetFakeRatio(track->GetFakeRatio());
3338 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3339 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3340 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3344 delete forwardtrack;
3346 for (Int_t i=0;i<entries;i++){
3347 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3349 if (!track) continue;
3351 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3352 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3353 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3354 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3355 delete array->RemoveAt(i);
3365 SortTrackHypothesys(esdindex,checkmax,1);
3367 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3368 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3369 besttrack = (AliITStrackMI*)array->At(0);
3370 if (!besttrack) return 0;
3371 besttrack->SetChi2MIP(8,0);
3372 fBestTrackIndex[esdindex]=0;
3373 entries = array->GetEntriesFast();
3374 AliITStrackMI *longtrack =0;
3376 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3377 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3378 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3379 if (!track->GetConstrain()) continue;
3380 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3381 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3382 if (track->GetChi2MIP(0)>4.) continue;
3383 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3386 //if (longtrack) besttrack=longtrack;
3389 AliITSRecPoint * clist[6];
3390 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3391 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3392 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3393 RegisterClusterTracks(besttrack,esdindex);
3400 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3401 if (sharedtrack>=0){
3403 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3405 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3411 if (shared>2.5) return 0;
3412 if (shared>1.0) return besttrack;
3414 // Don't sign clusters if not gold track
3416 if (!besttrack->IsGoldPrimary()) return besttrack;
3417 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3419 if (fConstraint[fPass]){
3423 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3424 for (Int_t i=0;i<6;i++){
3425 Int_t index = besttrack->GetClIndex(i);
3426 if (index<=0) continue;
3427 Int_t ilayer = (index & 0xf0000000) >> 28;
3428 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3429 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3431 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3432 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3433 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3434 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3435 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3436 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3438 Bool_t cansign = kTRUE;
3439 for (Int_t itrack=0;itrack<entries; itrack++){
3440 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3441 if (!track) continue;
3442 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3443 if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3449 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3450 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3451 if (!c->IsUsed()) c->Use();
3457 //------------------------------------------------------------------------
3458 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3461 // get "best" hypothesys
3464 Int_t nentries = itsTracks.GetEntriesFast();
3465 for (Int_t i=0;i<nentries;i++){
3466 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3467 if (!track) continue;
3468 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3469 if (!array) continue;
3470 if (array->GetEntriesFast()<=0) continue;
3472 AliITStrackMI* longtrack=0;
3474 Float_t maxchi2=1000;
3475 for (Int_t j=0;j<array->GetEntriesFast();j++){
3476 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3477 if (!trackHyp) continue;
3478 if (trackHyp->GetGoldV0()) {
3479 longtrack = trackHyp; //gold V0 track taken
3482 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3483 Float_t chi2 = trackHyp->GetChi2MIP(0);
3485 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3487 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3489 if (chi2 > maxchi2) continue;
3490 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3497 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3498 if (!longtrack) {longtrack = besttrack;}
3499 else besttrack= longtrack;
3503 AliITSRecPoint * clist[6];
3504 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3506 track->SetNUsed(shared);
3507 track->SetNSkipped(besttrack->GetNSkipped());
3508 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3510 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3514 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3515 //if (sharedtrack==-1) sharedtrack=0;
3516 if (sharedtrack>=0) {
3517 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3520 if (besttrack&&fAfterV0) {
3521 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3523 if (besttrack&&fConstraint[fPass])
3524 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3525 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3526 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3527 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3533 //------------------------------------------------------------------------
3534 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3535 //--------------------------------------------------------------------
3536 //This function "cooks" a track label. If label<0, this track is fake.
3537 //--------------------------------------------------------------------
3540 if ( track->GetESDtrack()) tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3542 track->SetChi2MIP(9,0);
3544 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3545 Int_t cindex = track->GetClusterIndex(i);
3546 Int_t l=(cindex & 0xf0000000) >> 28;
3547 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3549 for (Int_t ind=0;ind<3;ind++){
3551 if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3552 AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3554 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3557 Int_t nclusters = track->GetNumberOfClusters();
3558 if (nclusters > 0) //PH Some tracks don't have any cluster
3559 track->SetFakeRatio(double(nwrong)/double(nclusters));
3561 if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3563 track->SetLabel(tpcLabel);
3565 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3568 //------------------------------------------------------------------------
3569 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3574 //AliITSRecPoint * clist[6];
3575 // Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3578 track->SetChi2MIP(9,0);
3579 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3580 Int_t cindex = track->GetClusterIndex(i);
3581 Int_t l=(cindex & 0xf0000000) >> 28;
3582 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3583 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3585 for (Int_t ind=0;ind<3;ind++){
3586 if (cl->GetLabel(ind)==lab) isWrong=0;
3588 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3590 //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue; //shared track
3591 //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3592 //if (l<4&& !(cl->GetType()==1)) continue;
3593 dedx[accepted]= track->GetSampledEdx(i);
3594 //dedx[accepted]= track->fNormQ[l];
3602 TMath::Sort(accepted,dedx,indexes,kFALSE);
3605 Double_t nl=low*accepted, nu =up*accepted;
3607 Float_t sumweight =0;
3608 for (Int_t i=0; i<accepted; i++) {
3610 if (i<nl+0.1) weight = TMath::Max(1.-(nl-i),0.);
3611 if (i>nu-1) weight = TMath::Max(nu-i,0.);
3612 sumamp+= dedx[indexes[i]]*weight;
3615 track->SetdEdx(sumamp/sumweight);
3617 //------------------------------------------------------------------------
3618 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3621 if (fCoefficients) delete []fCoefficients;
3622 fCoefficients = new Float_t[ntracks*48];
3623 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3625 //------------------------------------------------------------------------
3626 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3632 Float_t theta = track->GetTgl();
3633 Float_t phi = track->GetSnp();
3634 phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
3635 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3636 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3637 // Take into account the mis-alignment (bring track to cluster plane)
3638 Double_t xTrOrig=track->GetX();
3639 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3640 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3641 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3642 // Bring the track back to detector plane in ideal geometry
3643 // [mis-alignment will be accounted for in UpdateMI()]
3644 if (!track->Propagate(xTrOrig)) return 1000.;
3646 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3647 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3649 chi2+=0.5*TMath::Min(delta/2,2.);
3650 chi2+=2.*cluster->GetDeltaProbability();
3653 track->SetNy(layer,ny);
3654 track->SetNz(layer,nz);
3655 track->SetSigmaY(layer,erry);
3656 track->SetSigmaZ(layer, errz);
3657 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3658 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/(1.- track->GetSnp()*track->GetSnp())));
3662 //------------------------------------------------------------------------
3663 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3668 Int_t layer = (index & 0xf0000000) >> 28;
3669 track->SetClIndex(layer, index);
3670 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3671 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3672 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3673 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3677 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
3680 // Take into account the mis-alignment (bring track to cluster plane)
3681 Double_t xTrOrig=track->GetX();
3682 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3683 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3684 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3685 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3687 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3691 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3692 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3695 Int_t updated = track->UpdateMI(&c,chi2,index);
3697 // Bring the track back to detector plane in ideal geometry
3698 if (!track->Propagate(xTrOrig)) return 0;
3700 if(!updated) AliDebug(2,"update failed");
3704 //------------------------------------------------------------------------
3705 void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3708 //DCA sigmas parameterization
3709 //to be paramterized using external parameters in future
3712 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3713 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3715 //------------------------------------------------------------------------
3716 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3720 Int_t entries = ClusterArray->GetEntriesFast();
3721 if (entries<4) return;
3722 AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3723 Int_t layer = cluster->GetLayer();
3724 if (layer>1) return;
3726 Int_t ncandidates=0;
3727 Float_t r = (layer>0)? 7:4;
3729 for (Int_t i=0;i<entries;i++){
3730 AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3731 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3732 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3733 index[ncandidates] = i; //candidate to belong to delta electron track
3735 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3736 cl0->SetDeltaProbability(1);
3742 for (Int_t i=0;i<ncandidates;i++){
3743 AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3744 if (cl0->GetDeltaProbability()>0.8) continue;
3747 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3748 sumy=sumz=sumy2=sumyz=sumw=0.0;
3749 for (Int_t j=0;j<ncandidates;j++){
3751 AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3753 Float_t dz = cl0->GetZ()-cl1->GetZ();
3754 Float_t dy = cl0->GetY()-cl1->GetY();
3755 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3756 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3757 y[ncl] = cl1->GetY();
3758 z[ncl] = cl1->GetZ();
3759 sumy+= y[ncl]*weight;
3760 sumz+= z[ncl]*weight;
3761 sumy2+=y[ncl]*y[ncl]*weight;
3762 sumyz+=y[ncl]*z[ncl]*weight;
3767 if (ncl<4) continue;
3768 Float_t det = sumw*sumy2 - sumy*sumy;
3770 if (TMath::Abs(det)>0.01){
3771 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3772 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3773 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3776 Float_t z0 = sumyz/sumy;
3777 delta = TMath::Abs(cl0->GetZ()-z0);
3780 cl0->SetDeltaProbability(1-20.*delta);
3784 //------------------------------------------------------------------------
3785 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3789 track->UpdateESDtrack(flags);
3790 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3791 if (oldtrack) delete oldtrack;
3792 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3793 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3794 printf("Problem\n");
3797 //------------------------------------------------------------------------
3798 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3800 // Get nearest upper layer close to the point xr.
3801 // rough approximation
3803 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3804 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3806 for (Int_t i=0;i<6;i++){
3807 if (radius<kRadiuses[i]){
3814 //------------------------------------------------------------------------
3815 void AliITStrackerMI::UpdateTPCV0(AliESDEvent *event){
3817 //try to update, or reject TPC V0s
3819 Int_t nv0s = event->GetNumberOfV0s();
3820 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3822 for (Int_t i=0;i<nv0s;i++){
3823 AliESDv0 * vertex = event->GetV0(i);
3824 Int_t ip = vertex->GetIndex(0);
3825 Int_t im = vertex->GetIndex(1);
3827 TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3828 TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3829 AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3830 AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3834 if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3835 if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3836 if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3841 if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3842 if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3843 if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3846 if (vertex->GetStatus()==-100) continue;
3848 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
3849 Int_t clayer = GetNearestLayer(xrp); //I.B.
3850 vertex->SetNBefore(clayer); //
3851 vertex->SetChi2Before(9*clayer); //
3852 vertex->SetNAfter(6-clayer); //
3853 vertex->SetChi2After(0); //
3855 if (clayer >1 ){ // calculate chi2 before vertex
3856 Float_t chi2p = 0, chi2m=0;
3859 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3860 if (trackp->GetClIndex(ilayer)>0){
3861 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3862 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3873 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3874 if (trackm->GetClIndex(ilayer)>0){
3875 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3876 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3885 vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3886 if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex
3889 if (clayer < 5 ){ // calculate chi2 after vertex
3890 Float_t chi2p = 0, chi2m=0;
3892 if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3893 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3894 if (trackp->GetClIndex(ilayer)>0){
3895 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3896 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3906 if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3907 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3908 if (trackm->GetClIndex(ilayer)>0){
3909 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3910 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3919 vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3920 if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS
3925 //------------------------------------------------------------------------
3926 void AliITStrackerMI::FindV02(AliESDEvent *event)
3931 // Cuts on DCA - R dependent
3932 // max distance DCA between 2 tracks cut
3933 // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3935 const Float_t kMaxDist0 = 0.1;
3936 const Float_t kMaxDist1 = 0.1;
3937 const Float_t kMaxDist = 1;
3938 const Float_t kMinPointAngle = 0.85;
3939 const Float_t kMinPointAngle2 = 0.99;
3940 const Float_t kMinR = 0.5;
3941 const Float_t kMaxR = 220;
3942 //const Float_t kCausality0Cut = 0.19;
3943 //const Float_t kLikelihood01Cut = 0.25;
3944 //const Float_t kPointAngleCut = 0.9996;
3945 const Float_t kCausality0Cut = 0.19;
3946 const Float_t kLikelihood01Cut = 0.45;
3947 const Float_t kLikelihood1Cut = 0.5;
3948 const Float_t kCombinedCut = 0.55;
3952 TTreeSRedirector &cstream = *fDebugStreamer;
3953 Int_t ntracks = event->GetNumberOfTracks();
3954 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3955 fOriginal.Expand(ntracks);
3956 fTrackHypothesys.Expand(ntracks);
3957 fBestHypothesys.Expand(ntracks);
3959 AliHelix * helixes = new AliHelix[ntracks+2];
3960 TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain
3961 TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain
3962 TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain
3963 Bool_t * forbidden = new Bool_t [ntracks+2];
3964 Int_t *itsmap = new Int_t [ntracks+2];
3965 Float_t *dist = new Float_t[ntracks+2];
3966 Float_t *normdist0 = new Float_t[ntracks+2];
3967 Float_t *normdist1 = new Float_t[ntracks+2];
3968 Float_t *normdist = new Float_t[ntracks+2];
3969 Float_t *norm = new Float_t[ntracks+2];
3970 Float_t *maxr = new Float_t[ntracks+2];
3971 Float_t *minr = new Float_t[ntracks+2];
3972 Float_t *minPointAngle= new Float_t[ntracks+2];
3974 AliV0 *pvertex = new AliV0;
3975 AliITStrackMI * dummy= new AliITStrackMI;
3977 AliITStrackMI trackat0; //temporary track for DCA calculation
3979 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3981 // make ITS - ESD map
3983 for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3984 itsmap[itrack] = -1;
3985 forbidden[itrack] = kFALSE;
3986 maxr[itrack] = kMaxR;
3987 minr[itrack] = kMinR;
3988 minPointAngle[itrack] = kMinPointAngle;
3990 for (Int_t itrack=0;itrack<nitstracks;itrack++){
3991 AliITStrackMI * original = (AliITStrackMI*)(fOriginal.At(itrack));
3992 Int_t esdindex = original->GetESDtrack()->GetID();
3993 itsmap[esdindex] = itrack;
3996 // create ITS tracks from ESD tracks if not done before
3998 for (Int_t itrack=0;itrack<ntracks;itrack++){
3999 if (itsmap[itrack]>=0) continue;
4000 AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
4001 //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
4002 //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
4003 tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP()); //I.B.
4004 if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
4005 // tracks which can reach inner part of ITS
4006 // propagate track to outer its volume - with correction for material
4007 CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
4009 itsmap[itrack] = nitstracks;
4010 fOriginal.AddAt(tpctrack,nitstracks);
4014 // fill temporary arrays
4016 for (Int_t itrack=0;itrack<ntracks;itrack++){
4017 AliESDtrack * esdtrack = event->GetTrack(itrack);
4018 Int_t itsindex = itsmap[itrack];
4019 AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
4020 if (!original) continue;
4021 AliITStrackMI *bestConst = 0;
4022 AliITStrackMI *bestLong = 0;
4023 AliITStrackMI *best = 0;
4026 TObjArray * array = (TObjArray*) fTrackHypothesys.At(itsindex);
4027 Int_t hentries = (array==0) ? 0 : array->GetEntriesFast();
4028 // Get best track with vertex constrain
4029 for (Int_t ih=0;ih<hentries;ih++){
4030 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4031 if (!trackh->GetConstrain()) continue;
4032 if (!bestConst) bestConst = trackh;
4033 if (trackh->GetNumberOfClusters()>5.0){
4034 bestConst = trackh; // full track - with minimal chi2
4037 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()) continue;
4041 // Get best long track without vertex constrain and best track without vertex constrain
4042 for (Int_t ih=0;ih<hentries;ih++){
4043 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4044 if (trackh->GetConstrain()) continue;
4045 if (!best) best = trackh;
4046 if (!bestLong) bestLong = trackh;
4047 if (trackh->GetNumberOfClusters()>5.0){
4048 bestLong = trackh; // full track - with minimal chi2
4051 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()) continue;
4056 bestLong = original;
4058 //I.B. trackat0 = *bestLong;
4059 new (&trackat0) AliITStrackMI(*bestLong);
4060 Double_t xx,yy,zz,alpha;
4061 if (!bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz)) continue;
4062 alpha = TMath::ATan2(yy,xx);
4063 if (!trackat0.Propagate(alpha,0)) continue;
4064 // calculate normalized distances to the vertex
4066 Float_t ptfac = (1.+100.*TMath::Abs(trackat0.GetC()));
4067 if ( bestLong->GetNumberOfClusters()>3 ){
4068 dist[itsindex] = trackat0.GetY();
4069 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4070 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4071 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4072 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4074 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
4075 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
4076 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
4078 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
4079 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
4083 if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
4084 dist[itsindex] = bestConst->GetD(0);
4085 norm[itsindex] = bestConst->GetDnorm(0);
4086 normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4087 normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4088 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4090 dist[itsindex] = trackat0.GetY();
4091 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4092 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4093 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4094 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4095 if (TMath::Abs(trackat0.GetTgl())>1.05){
4096 if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
4097 if (normdist[itsindex]>3) {
4098 minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
4104 //-----------------------------------------------------------
4105 //Forbid primary track candidates -
4107 //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
4108 //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
4109 //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
4110 //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
4111 //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
4112 //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
4113 //-----------------------------------------------------------
4115 if (bestLong->GetNumberOfClusters()<4 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5) forbidden[itsindex]=kTRUE;
4116 if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5) forbidden[itsindex]=kTRUE;
4117 if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
4118 if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0) forbidden[itsindex]=kTRUE;
4119 if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2) forbidden[itsindex]=kTRUE;
4120 if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1) forbidden[itsindex]=kTRUE;
4121 if (bestConst->GetNormChi2(0)<2.5) {
4122 minPointAngle[itsindex]= 0.9999;
4123 maxr[itsindex] = 10;
4127 //forbid daughter kink candidates
4129 if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
4130 Bool_t isElectron = kTRUE;
4131 Bool_t isProton = kTRUE;
4133 esdtrack->GetESDpid(pid);
4134 for (Int_t i=1;i<5;i++){
4135 if (pid[0]<pid[i]) isElectron= kFALSE;
4136 if (pid[4]<pid[i]) isProton= kFALSE;
4139 forbidden[itsindex]=kFALSE;
4140 normdist[itsindex]*=-1;
4143 if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
4144 normdist[itsindex]*=-1;
4148 // Causality cuts in TPC volume
4150 if (esdtrack->GetTPCdensity(0,10) >0.6) maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
4151 if (esdtrack->GetTPCdensity(10,30)>0.6) maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
4152 if (esdtrack->GetTPCdensity(20,40)>0.6) maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
4153 if (esdtrack->GetTPCdensity(30,50)>0.6) maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
4155 if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
4161 "Tr1.="<<((bestConst)? bestConst:dummy)<<
4163 "Tr3.="<<&trackat0<<
4165 "Dist="<<dist[itsindex]<<
4166 "ND0="<<normdist0[itsindex]<<
4167 "ND1="<<normdist1[itsindex]<<
4168 "ND="<<normdist[itsindex]<<
4169 "Pz="<<primvertex[2]<<
4170 "Forbid="<<forbidden[itsindex]<<
4174 trackarray.AddAt(best,itsindex);
4175 trackarrayc.AddAt(bestConst,itsindex);
4176 trackarrayl.AddAt(bestLong,itsindex);
4177 new (&helixes[itsindex]) AliHelix(*best);
4182 // first iterration of V0 finder
4184 for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
4185 Int_t itrack0 = itsmap[iesd0];
4186 if (forbidden[itrack0]) continue;
4187 AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
4188 if (!btrack0) continue;
4189 if (btrack0->GetSign()>0) continue;
4190 AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
4192 for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
4193 Int_t itrack1 = itsmap[iesd1];
4194 if (forbidden[itrack1]) continue;
4196 AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
4197 if (!btrack1) continue;
4198 if (btrack1->GetSign()<0) continue;
4199 Bool_t isGold = kFALSE;
4200 if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
4203 AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
4204 AliHelix &h1 = helixes[itrack0];
4205 AliHelix &h2 = helixes[itrack1];
4207 // find linear distance
4212 Double_t phase[2][2],radius[2];
4213 Int_t points = h1.GetRPHIintersections(h2, phase, radius);
4214 if (points==0) continue;
4215 Double_t delta[2]={1000000,1000000};
4217 h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
4219 if (radius[1]<rmin) rmin = radius[1];
4220 h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
4222 rmin = TMath::Sqrt(rmin);
4223 Double_t distance = 0;
4224 Double_t radiusC = 0;
4226 if (points==1 || delta[0]<delta[1]){
4227 distance = TMath::Sqrt(delta[0]);
4228 radiusC = TMath::Sqrt(radius[0]);
4230 distance = TMath::Sqrt(delta[1]);
4231 radiusC = TMath::Sqrt(radius[1]);
4234 if (radiusC<TMath::Max(minr[itrack0],minr[itrack1])) continue;
4235 if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1])) continue;
4236 Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
4237 if (distance>maxDist) continue;
4238 Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
4239 if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
4242 // Double_t distance = TestV0(h1,h2,pvertex,rmin);
4244 // if (distance>maxDist) continue;
4245 // if (pvertex->GetRr()<kMinR) continue;
4246 // if (pvertex->GetRr()>kMaxR) continue;
4247 AliITStrackMI * track0=btrack0;
4248 AliITStrackMI * track1=btrack1;
4249 // if (pvertex->GetRr()<3.5){
4251 //use longest tracks inside the pipe
4252 track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4253 track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4257 pvertex->SetParamN(*track0);
4258 pvertex->SetParamP(*track1);
4259 pvertex->Update(primvertex);
4260 pvertex->SetClusters(track0->ClIndex(),track1->ClIndex()); // register clusters
4262 if (pvertex->GetRr()<kMinR) continue;
4263 if (pvertex->GetRr()>kMaxR) continue;
4264 if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4265 //Bo: if (pvertex->GetDist2()>maxDist) continue;
4266 if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4267 //Bo: pvertex->SetLab(0,track0->GetLabel());
4268 //Bo: pvertex->SetLab(1,track1->GetLabel());
4269 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4270 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4272 AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4273 AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4277 TObjArray * array0b = (TObjArray*)fBestHypothesys.At(itrack0);
4278 if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4279 fCurrentEsdTrack = itrack0;
4280 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4282 TObjArray * array1b = (TObjArray*)fBestHypothesys.At(itrack1);
4283 if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4284 fCurrentEsdTrack = itrack1;
4285 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4288 AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4289 AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4290 AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4291 AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4293 Float_t minchi2before0=16;
4294 Float_t minchi2before1=16;
4295 Float_t minchi2after0 =16;
4296 Float_t minchi2after1 =16;
4297 Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
4298 Int_t maxLayer = GetNearestLayer(xrp); //I.B.
4300 if (array0b) for (Int_t i=0;i<5;i++){
4301 // best track after vertex
4302 AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4303 if (!btrack) continue;
4304 if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4305 // if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4306 if (btrack->GetX()<pvertex->GetRr()-2.) {
4307 if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4310 if (maxLayer<3){ // take prim vertex as additional measurement
4311 if (normdist[itrack0]>0 && htrackc0){
4312 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4314 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4318 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4320 if (!btrack->GetClIndex(ilayer)){
4324 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4325 for (Int_t itrack=0;itrack<4;itrack++){
4326 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4327 sumchi2+=18.; //shared cluster
4331 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4332 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4336 if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4338 continue; //safety space - Geo manager will give exact layer
4341 minchi2after0 = btrack->GetNormChi2(i);
4344 if (array1b) for (Int_t i=0;i<5;i++){
4345 // best track after vertex
4346 AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4347 if (!btrack) continue;
4348 if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4349 // if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4350 if (btrack->GetX()<pvertex->GetRr()-2){
4351 if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4354 if (maxLayer<3){ // take prim vertex as additional measurement
4355 if (normdist[itrack1]>0 && htrackc1){
4356 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4358 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4362 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4364 if (!btrack->GetClIndex(ilayer)){
4368 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4369 for (Int_t itrack=0;itrack<4;itrack++){
4370 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4371 sumchi2+=18.; //shared cluster
4375 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4376 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4380 if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4382 continue; //safety space - Geo manager will give exact layer
4385 minchi2after1 = btrack->GetNormChi2(i);
4389 // position resolution - used for DCA cut
4390 Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4391 (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4392 (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4393 sigmad =TMath::Sqrt(sigmad)+0.04;
4394 if (pvertex->GetRr()>50){
4395 Double_t cov0[15],cov1[15];
4396 track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4397 track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4398 sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4399 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4400 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4401 sigmad =TMath::Sqrt(sigmad)+0.05;
4405 vertex2.SetParamN(*track0b);
4406 vertex2.SetParamP(*track1b);
4407 vertex2.Update(primvertex);
4408 //Bo: if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4409 if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4410 pvertex->SetParamN(*track0b);
4411 pvertex->SetParamP(*track1b);
4412 pvertex->Update(primvertex);
4413 pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex()); // register clusters
4414 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4415 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4417 pvertex->SetDistSigma(sigmad);
4418 //Bo: pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4419 pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4421 // define likelihhod and causalities
4423 Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4425 Float_t fnorm0 = normdist[itrack0];
4426 if (fnorm0<0) fnorm0*=-3;
4427 Float_t fnorm1 = normdist[itrack1];
4428 if (fnorm1<0) fnorm1*=-3;
4429 if ((pvertex->GetAnglep()[2]>0.1) || ( (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 ) || (pvertex->GetRr()<3)){
4430 pb0 = TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4431 pb1 = TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4433 pvertex->SetChi2Before(normdist[itrack0]);
4434 pvertex->SetChi2After(normdist[itrack1]);
4435 pvertex->SetNAfter(0);
4436 pvertex->SetNBefore(0);
4438 pvertex->SetChi2Before(minchi2before0);
4439 pvertex->SetChi2After(minchi2before1);
4440 if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4441 pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4442 pb1 = TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4444 pvertex->SetNAfter(maxLayer);
4445 pvertex->SetNBefore(maxLayer);
4447 if (pvertex->GetRr()<90){
4448 pa0 *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4449 pa1 *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4451 if (pvertex->GetRr()<20){
4452 pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4453 pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4456 pvertex->SetCausality(pb0,pb1,pa0,pa1);
4458 // Likelihood calculations - apply cuts
4460 Bool_t v0OK = kTRUE;
4461 Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4462 p12 += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4463 p12 = TMath::Sqrt(p12); // "mean" momenta
4464 Float_t sigmap0 = 0.0001+0.001/(0.1+pvertex->GetRr());
4465 Float_t sigmap = 0.5*sigmap0*(0.6+0.4*p12); // "resolution: of point angle - as a function of radius and momenta
4467 Float_t causalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4468 Float_t causalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4469 TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4471 //Bo: Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4472 Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4473 Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4475 Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4476 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4477 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4478 0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4480 if (causalityA<kCausality0Cut) v0OK = kFALSE;
4481 if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut) v0OK = kFALSE;
4482 if (likelihood1<kLikelihood1Cut) v0OK = kFALSE;
4483 if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4488 Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4490 "Tr0.="<<track0<< //best without constrain
4491 "Tr1.="<<track1<< //best without constrain
4492 "Tr0B.="<<track0b<< //best without constrain after vertex
4493 "Tr1B.="<<track1b<< //best without constrain after vertex
4494 "Tr0C.="<<htrackc0<< //best with constrain if exist
4495 "Tr1C.="<<htrackc1<< //best with constrain if exist
4496 "Tr0L.="<<track0l<< //longest best
4497 "Tr1L.="<<track1l<< //longest best
4498 "Esd0.="<<track0->GetESDtrack()<< // esd track0 params
4499 "Esd1.="<<track1->GetESDtrack()<< // esd track1 params
4500 "V0.="<<pvertex<< //vertex properties
4501 "V0b.="<<&vertex2<< //vertex properties at "best" track
4502 "ND0="<<normdist[itrack0]<< //normalize distance for track0
4503 "ND1="<<normdist[itrack1]<< //normalize distance for track1
4505 // "RejectBase="<<rejectBase<< //rejection in First itteration
4511 //if (rejectBase) continue;
4513 pvertex->SetStatus(0);
4514 // if (rejectBase) {
4515 // pvertex->SetStatus(-100);
4517 if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4518 //Bo: pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4519 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4520 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4522 // AliV0vertex vertexjuri(*track0,*track1);
4523 // vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4524 // event->AddV0(&vertexjuri);
4525 pvertex->SetStatus(100);
4527 pvertex->SetOnFlyStatus(kTRUE);
4528 pvertex->ChangeMassHypothesis(kK0Short);
4529 event->AddV0(pvertex);
4535 // delete temporary arrays
4538 delete[] minPointAngle;
4550 //------------------------------------------------------------------------
4551 void AliITStrackerMI::RefitV02(AliESDEvent *event)
4554 //try to refit V0s in the third path of the reconstruction
4556 TTreeSRedirector &cstream = *fDebugStreamer;
4558 Int_t nv0s = event->GetNumberOfV0s();
4559 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4561 for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4562 AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4563 if (!v0mi) continue;
4564 Int_t itrack0 = v0mi->GetIndex(0);
4565 Int_t itrack1 = v0mi->GetIndex(1);
4566 AliESDtrack *esd0 = event->GetTrack(itrack0);
4567 AliESDtrack *esd1 = event->GetTrack(itrack1);
4568 if (!esd0||!esd1) continue;
4569 AliITStrackMI tpc0(*esd0);
4570 AliITStrackMI tpc1(*esd1);
4571 Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4572 Double_t alpha =TMath::ATan2(y,x); //I.B.
4573 if (v0mi->GetRr()>85){
4574 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4575 v0temp.SetParamN(tpc0);
4576 v0temp.SetParamP(tpc1);
4577 v0temp.Update(primvertex);
4578 if (kFALSE) cstream<<"Refit"<<
4580 "V0refit.="<<&v0temp<<
4584 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4585 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4586 v0mi->SetParamN(tpc0);
4587 v0mi->SetParamP(tpc1);
4588 v0mi->Update(primvertex);
4593 if (v0mi->GetRr()>35){
4594 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4595 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4596 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4597 v0temp.SetParamN(tpc0);
4598 v0temp.SetParamP(tpc1);
4599 v0temp.Update(primvertex);
4600 if (kFALSE) cstream<<"Refit"<<
4602 "V0refit.="<<&v0temp<<
4606 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4607 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4608 v0mi->SetParamN(tpc0);
4609 v0mi->SetParamP(tpc1);
4610 v0mi->Update(primvertex);
4615 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4616 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4617 // if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4618 if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4619 v0temp.SetParamN(tpc0);
4620 v0temp.SetParamP(tpc1);
4621 v0temp.Update(primvertex);
4622 if (kFALSE) cstream<<"Refit"<<
4624 "V0refit.="<<&v0temp<<
4628 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4629 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4630 v0mi->SetParamN(tpc0);
4631 v0mi->SetParamP(tpc1);
4632 v0mi->Update(primvertex);
4637 //------------------------------------------------------------------------
4638 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4639 //--------------------------------------------------------------------
4640 // Fill a look-up table with mean material
4641 //--------------------------------------------------------------------
4645 Double_t point1[3],point2[3];
4646 Double_t phi,cosphi,sinphi,z;
4647 // 0-5 layers, 6 pipe, 7-8 shields
4648 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4649 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4651 Int_t ifirst=0,ilast=0;
4652 if(material.Contains("Pipe")) {
4654 } else if(material.Contains("Shields")) {
4656 } else if(material.Contains("Layers")) {
4659 Error("BuildMaterialLUT","Wrong layer name\n");
4662 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4663 Double_t param[5]={0.,0.,0.,0.,0.};
4664 for (Int_t i=0; i<n; i++) {
4665 phi = 2.*TMath::Pi()*gRandom->Rndm();
4666 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4667 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4668 point1[0] = rmin[imat]*cosphi;
4669 point1[1] = rmin[imat]*sinphi;
4671 point2[0] = rmax[imat]*cosphi;
4672 point2[1] = rmax[imat]*sinphi;
4674 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4675 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4677 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4679 fxOverX0Layer[imat] = param[1];
4680 fxTimesRhoLayer[imat] = param[0]*param[4];
4681 } else if(imat==6) {
4682 fxOverX0Pipe = param[1];
4683 fxTimesRhoPipe = param[0]*param[4];
4684 } else if(imat==7) {
4685 fxOverX0Shield[0] = param[1];
4686 fxTimesRhoShield[0] = param[0]*param[4];
4687 } else if(imat==8) {
4688 fxOverX0Shield[1] = param[1];
4689 fxTimesRhoShield[1] = param[0]*param[4];
4693 printf("%s\n",material.Data());
4694 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4695 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4696 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4697 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4698 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4699 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4700 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4701 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4702 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4706 //------------------------------------------------------------------------
4707 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4708 TString direction) {
4709 //-------------------------------------------------------------------
4710 // Propagate beyond beam pipe and correct for material
4711 // (material budget in different ways according to fUseTGeo value)
4712 //-------------------------------------------------------------------
4714 // Define budget mode:
4715 // 0: material from AliITSRecoParam (hard coded)
4716 // 1: material from TGeo in one step (on the fly)
4717 // 2: material from lut
4718 // 3: material from TGeo in one step (same for all hypotheses)
4731 if(fTrackingPhase.Contains("Clusters2Tracks"))
4732 { mode=3; } else { mode=1; }
4735 if(fTrackingPhase.Contains("Clusters2Tracks"))
4736 { mode=3; } else { mode=2; }
4742 if(fTrackingPhase.Contains("Default")) mode=0;
4744 Int_t index=fCurrentEsdTrack;
4746 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4747 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4749 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4751 Double_t xOverX0,x0,lengthTimesMeanDensity;
4752 Bool_t anglecorr=kTRUE;
4756 xOverX0 = AliITSRecoParam::GetdPipe();
4757 x0 = AliITSRecoParam::GetX0Be();
4758 lengthTimesMeanDensity = xOverX0*x0;
4759 lengthTimesMeanDensity *= dir;
4760 if (!t->Propagate(xToGo)) return 0;
4761 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4764 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4767 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4768 xOverX0 = fxOverX0Pipe;
4769 lengthTimesMeanDensity = fxTimesRhoPipe;
4770 lengthTimesMeanDensity *= dir;
4771 if (!t->Propagate(xToGo)) return 0;
4772 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4775 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4776 if(fxOverX0PipeTrks[index]<0) {
4777 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4778 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4779 (1.-t->GetSnp()*t->GetSnp()));
4780 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4781 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4784 xOverX0 = fxOverX0PipeTrks[index];
4785 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4786 lengthTimesMeanDensity *= dir;
4787 if (!t->Propagate(xToGo)) return 0;
4788 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4794 //------------------------------------------------------------------------
4795 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4797 TString direction) {
4798 //-------------------------------------------------------------------
4799 // Propagate beyond SPD or SDD shield and correct for material
4800 // (material budget in different ways according to fUseTGeo value)
4801 //-------------------------------------------------------------------
4803 // Define budget mode:
4804 // 0: material from AliITSRecoParam (hard coded)
4805 // 1: material from TGeo in steps of X cm (on the fly)
4806 // X = AliITSRecoParam::GetStepSizeTGeo()
4807 // 2: material from lut
4808 // 3: material from TGeo in one step (same for all hypotheses)
4821 if(fTrackingPhase.Contains("Clusters2Tracks"))
4822 { mode=3; } else { mode=1; }
4825 if(fTrackingPhase.Contains("Clusters2Tracks"))
4826 { mode=3; } else { mode=2; }
4832 if(fTrackingPhase.Contains("Default")) mode=0;
4834 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4836 Int_t shieldindex=0;
4837 if (shield.Contains("SDD")) { // SDDouter
4838 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4840 } else if (shield.Contains("SPD")) { // SPDouter
4841 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4844 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4848 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4850 Int_t index=2*fCurrentEsdTrack+shieldindex;
4852 Double_t xOverX0,x0,lengthTimesMeanDensity;
4853 Bool_t anglecorr=kTRUE;
4858 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4859 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4860 lengthTimesMeanDensity = xOverX0*x0;
4861 lengthTimesMeanDensity *= dir;
4862 if (!t->Propagate(xToGo)) return 0;
4863 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4866 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4867 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4870 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4871 xOverX0 = fxOverX0Shield[shieldindex];
4872 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4873 lengthTimesMeanDensity *= dir;
4874 if (!t->Propagate(xToGo)) return 0;
4875 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4878 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4879 if(fxOverX0ShieldTrks[index]<0) {
4880 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4881 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4882 (1.-t->GetSnp()*t->GetSnp()));
4883 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4884 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4887 xOverX0 = fxOverX0ShieldTrks[index];
4888 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4889 lengthTimesMeanDensity *= dir;
4890 if (!t->Propagate(xToGo)) return 0;
4891 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4897 //------------------------------------------------------------------------
4898 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4900 Double_t oldGlobXYZ[3],
4901 TString direction) {
4902 //-------------------------------------------------------------------
4903 // Propagate beyond layer and correct for material
4904 // (material budget in different ways according to fUseTGeo value)
4905 //-------------------------------------------------------------------
4907 // Define budget mode:
4908 // 0: material from AliITSRecoParam (hard coded)
4909 // 1: material from TGeo in stepsof X cm (on the fly)
4910 // X = AliITSRecoParam::GetStepSizeTGeo()
4911 // 2: material from lut
4912 // 3: material from TGeo in one step (same for all hypotheses)
4925 if(fTrackingPhase.Contains("Clusters2Tracks"))
4926 { mode=3; } else { mode=1; }
4929 if(fTrackingPhase.Contains("Clusters2Tracks"))
4930 { mode=3; } else { mode=2; }
4936 if(fTrackingPhase.Contains("Default")) mode=0;
4938 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4940 Double_t r=fgLayers[layerindex].GetR();
4941 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4943 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4945 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4947 Int_t index=6*fCurrentEsdTrack+layerindex;
4950 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4952 Bool_t anglecorr=kTRUE;
4956 Bool_t addTime = kFALSE;
4959 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4960 lengthTimesMeanDensity = xOverX0*x0;
4961 // Bring the track beyond the material
4962 if (!t->Propagate(xToGo)) return 0;
4963 lengthTimesMeanDensity *= dir;
4964 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4967 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4968 if (!t->GetLocalXat(rOld,xOld)) return 0;
4969 if (!t->Propagate(xOld)) return 0; // back before material (no correction)
4970 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4971 if (!t->PropagateToTGeo(xToGo,nsteps,addTime)) return 0; // cross the material and apply correction
4974 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4975 xOverX0 = fxOverX0Layer[layerindex];
4976 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4977 // Bring the track beyond the material
4978 if (!t->Propagate(xToGo)) return 0;
4979 lengthTimesMeanDensity *= dir;
4980 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4983 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4984 // Bring the track beyond the material
4985 if (!t->Propagate(xToGo)) return 0;
4986 Double_t globXYZ[3];
4987 if (!t->GetXYZ(globXYZ)) return 0;
4988 if (fxOverX0LayerTrks[index]<0) {
4989 AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4990 if(mparam[1]>900000) return 0;
4991 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4992 (1.-t->GetSnp()*t->GetSnp()));
4993 xOverX0=mparam[1]/angle;
4994 lengthTimesMeanDensity=mparam[0]*mparam[4]/angle;
4995 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0);
4996 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity);
4998 xOverX0 = fxOverX0LayerTrks[index];
4999 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
5000 lengthTimesMeanDensity *= dir;
5001 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
5007 //------------------------------------------------------------------------
5008 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
5009 //-----------------------------------------------------------------
5010 // Initialize LUT for storing material for each prolonged track
5011 //-----------------------------------------------------------------
5012 fxOverX0PipeTrks = new Float_t[ntracks];
5013 fxTimesRhoPipeTrks = new Float_t[ntracks];
5014 fxOverX0ShieldTrks = new Float_t[ntracks*2];
5015 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
5016 fxOverX0LayerTrks = new Float_t[ntracks*6];
5017 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
5019 for(Int_t i=0; i<ntracks; i++) {
5020 fxOverX0PipeTrks[i] = -1.;
5021 fxTimesRhoPipeTrks[i] = -1.;
5023 for(Int_t j=0; j<ntracks*2; j++) {
5024 fxOverX0ShieldTrks[j] = -1.;
5025 fxTimesRhoShieldTrks[j] = -1.;
5027 for(Int_t k=0; k<ntracks*6; k++) {
5028 fxOverX0LayerTrks[k] = -1.;
5029 fxTimesRhoLayerTrks[k] = -1.;
5036 //------------------------------------------------------------------------
5037 void AliITStrackerMI::DeleteTrksMaterialLUT() {
5038 //-----------------------------------------------------------------
5039 // Delete LUT for storing material for each prolonged track
5040 //-----------------------------------------------------------------
5041 if(fxOverX0PipeTrks) {
5042 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
5044 if(fxOverX0ShieldTrks) {
5045 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
5048 if(fxOverX0LayerTrks) {
5049 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
5051 if(fxTimesRhoPipeTrks) {
5052 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
5054 if(fxTimesRhoShieldTrks) {
5055 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
5057 if(fxTimesRhoLayerTrks) {
5058 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
5062 //------------------------------------------------------------------------
5063 Int_t AliITStrackerMI::CheckSkipLayer(AliITStrackMI *track,
5064 Int_t ilayer,Int_t idet) const {
5065 //-----------------------------------------------------------------
5066 // This method is used to decide whether to allow a prolongation
5067 // without clusters, because we want to skip the layer.
5068 // In this case the return value is > 0:
5069 // return 1: the user requested to skip a layer
5070 // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
5071 //-----------------------------------------------------------------
5073 if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
5075 if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
5076 // check if track will cross SPD outer layer
5077 Double_t phiAtSPD2,zAtSPD2;
5078 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
5079 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
5085 //------------------------------------------------------------------------
5086 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
5087 Int_t ilayer,Int_t idet,
5088 Double_t dz,Double_t dy,
5089 Bool_t noClusters) const {
5090 //-----------------------------------------------------------------
5091 // This method is used to decide whether to allow a prolongation
5092 // without clusters, because there is a dead zone in the road.
5093 // In this case the return value is > 0:
5094 // return 1: dead zone at z=0,+-7cm in SPD
5095 // return 2: all road is "bad" (dead or noisy) from the OCDB
5096 // return 3: something "bad" (dead or noisy) from the OCDB
5097 //-----------------------------------------------------------------
5099 // check dead zones at z=0,+-7cm in the SPD
5100 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
5101 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5102 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5103 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
5104 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5105 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5106 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
5107 for (Int_t i=0; i<3; i++)
5108 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
5109 AliDebug(2,Form("crack SPD %d",ilayer));
5114 // check bad zones from OCDB
5115 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
5117 if (idet<0) return 0;
5119 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
5122 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
5123 if (ilayer==0 || ilayer==1) { // ---------- SPD
5125 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
5127 detSizeFactorX *= 2.;
5128 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
5131 AliITSsegmentation *segm = (AliITSsegmentation*)fDetTypeRec->GetSegmentationModel(detType);
5132 if (detType==2) segm->SetLayer(ilayer+1);
5133 Float_t detSizeX = detSizeFactorX*segm->Dx();
5134 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
5136 // check if the road overlaps with bad chips
5138 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
5139 Float_t zlocmin = zloc-dz;
5140 Float_t zlocmax = zloc+dz;
5141 Float_t xlocmin = xloc-dy;
5142 Float_t xlocmax = xloc+dy;
5143 Int_t chipsInRoad[100];
5145 // check if road goes out of detector
5146 Bool_t touchNeighbourDet=kFALSE;
5147 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
5148 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
5149 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5150 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
5151 AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
5153 // check if this detector is bad
5155 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
5156 if(!touchNeighbourDet) {
5157 return 2; // all detectors in road are bad
5159 return 3; // at least one is bad
5163 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
5164 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
5165 if (!nChipsInRoad) return 0;
5167 Bool_t anyBad=kFALSE,anyGood=kFALSE;
5168 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
5169 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
5170 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
5171 if (det.IsChipBad(chipsInRoad[iCh])) {
5179 if(!touchNeighbourDet) {
5180 AliDebug(2,"all bad in road");
5181 return 2; // all chips in road are bad
5183 return 3; // at least a bad chip in road
5188 AliDebug(2,"at least a bad in road");
5189 return 3; // at least a bad chip in road
5193 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
5194 || ilayer==4 || ilayer==5 // SSD
5195 || !noClusters) return 0;
5197 // There are no clusters in road: check if there is at least
5198 // a bad SPD pixel or SDD anode
5200 Int_t idetInITS=idet;
5201 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
5203 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
5204 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
5207 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
5211 //------------------------------------------------------------------------
5212 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
5213 AliITStrackMI *track,
5214 Float_t &xloc,Float_t &zloc) const {
5215 //-----------------------------------------------------------------
5216 // Gives position of track in local module ref. frame
5217 //-----------------------------------------------------------------
5222 if(idet<0) return kFALSE;
5224 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
5226 Int_t lad = Int_t(idet/ndet) + 1;
5228 Int_t det = idet - (lad-1)*ndet + 1;
5230 Double_t xyzGlob[3],xyzLoc[3];
5232 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
5233 // take into account the misalignment: xyz at real detector plane
5234 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
5236 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
5238 xloc = (Float_t)xyzLoc[0];
5239 zloc = (Float_t)xyzLoc[2];
5243 //------------------------------------------------------------------------
5244 Bool_t AliITStrackerMI::IsOKForPlaneEff(AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
5246 // Method to be optimized further:
5247 // Aim: decide whether a track can be used for PlaneEff evaluation
5248 // the decision is taken based on the track quality at the layer under study
5249 // no information on the clusters on this layer has to be used
5250 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
5251 // the cut is done on number of sigmas from the boundaries
5253 // Input: Actual track, layer [0,5] under study
5255 // Return: kTRUE if this is a good track
5257 // it will apply a pre-selection to obtain good quality tracks.
5258 // Here also you will have the possibility to put a control on the
5259 // impact point of the track on the basic block, in order to exclude border regions
5260 // this will be done by calling a proper method of the AliITSPlaneEff class.
5262 // input: AliITStrackMI* track, ilayer= layer number [0,5]
5263 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
5265 Int_t index[AliITSgeomTGeo::kNLayers];
5267 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
5269 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
5270 index[k]=clusters[k];
5274 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
5275 AliITSlayer &layer=fgLayers[ilayer];
5276 Double_t r=layer.GetR();
5277 AliITStrackMI tmp(*track);
5279 // require a minimal number of cluster in other layers and eventually clusters in closest layers
5281 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
5282 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
5283 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
5284 if (tmp.GetClIndex(lay)>0) ncl++;
5286 Bool_t nextout = kFALSE;
5287 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
5288 else nextout = ((tmp.GetClIndex(ilayer+1)>0)? kTRUE : kFALSE );
5289 Bool_t nextin = kFALSE;
5290 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
5291 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
5292 if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
5294 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
5295 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
5296 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
5297 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
5301 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
5302 Int_t idet=layer.FindDetectorIndex(phi,z);
5303 if(idet<0) { AliInfo(Form("cannot find detector"));
5306 // here check if it has good Chi Square.
5308 //propagate to the intersection with the detector plane
5309 const AliITSdetector &det=layer.GetDetector(idet);
5310 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
5314 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
5315 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5316 if(key>fPlaneEff->Nblock()) return kFALSE;
5317 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
5318 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
5320 // DEFINITION OF SEARCH ROAD FOR accepting a track
5322 //For the time being they are hard-wired, later on from AliITSRecoParam
5323 // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
5324 // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
5327 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5328 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5329 // done for RecPoints
5331 // exclude tracks at boundary between detectors
5332 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
5333 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
5334 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
5335 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
5336 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5338 if ( (locx-dx < blockXmn+boundaryWidth) ||
5339 (locx+dx > blockXmx-boundaryWidth) ||
5340 (locz-dz < blockZmn+boundaryWidth) ||
5341 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5344 //------------------------------------------------------------------------
5345 void AliITStrackerMI::UseTrackForPlaneEff(AliITStrackMI* track, Int_t ilayer) {
5347 // This Method has to be optimized! For the time-being it uses the same criteria
5348 // as those used in the search of extra clusters for overlapping modules.
5350 // Method Purpose: estabilish whether a track has produced a recpoint or not
5351 // in the layer under study (For Plane efficiency)
5353 // inputs: AliITStrackMI* track (pointer to a usable track)
5355 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5356 // traversed by this very track. In details:
5357 // - if a cluster can be associated to the track then call
5358 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5359 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5362 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5363 AliITSlayer &layer=fgLayers[ilayer];
5364 Double_t r=layer.GetR();
5365 AliITStrackMI tmp(*track);
5369 if (!tmp.GetPhiZat(r,phi,z)) return;
5370 Int_t idet=layer.FindDetectorIndex(phi,z);
5372 if(idet<0) { AliInfo(Form("cannot find detector"));
5376 //propagate to the intersection with the detector plane
5377 const AliITSdetector &det=layer.GetDetector(idet);
5378 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5382 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5384 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5385 TMath::Sqrt(tmp.GetSigmaZ2() +
5386 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5387 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5388 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5389 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5390 TMath::Sqrt(tmp.GetSigmaY2() +
5391 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5392 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5393 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5395 // road in global (rphi,z) [i.e. in tracking ref. system]
5396 Double_t zmin = tmp.GetZ() - dz;
5397 Double_t zmax = tmp.GetZ() + dz;
5398 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5399 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5401 // select clusters in road
5402 layer.SelectClusters(zmin,zmax,ymin,ymax);
5404 // Define criteria for track-cluster association
5405 Double_t msz = tmp.GetSigmaZ2() +
5406 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5407 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5408 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5409 Double_t msy = tmp.GetSigmaY2() +
5410 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5411 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5412 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5413 if (tmp.GetConstrain()) {
5414 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5415 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5417 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5418 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5420 msz = 1./msz; // 1/RoadZ^2
5421 msy = 1./msy; // 1/RoadY^2
5424 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5426 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5427 //Double_t tolerance=0.2;
5428 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5429 idetc = cl->GetDetectorIndex();
5430 if(idet!=idetc) continue;
5431 //Int_t ilay = cl->GetLayer();
5433 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5434 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5436 Double_t chi2=tmp.GetPredictedChi2(cl);
5437 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5441 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5443 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5444 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5445 if(key>fPlaneEff->Nblock()) return;
5446 Bool_t found=kFALSE;
5449 while ((cl=layer.GetNextCluster(clidx))!=0) {
5450 idetc = cl->GetDetectorIndex();
5451 if(idet!=idetc) continue;
5452 // here real control to see whether the cluster can be associated to the track.
5453 // cluster not associated to track
5454 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5455 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5456 // calculate track-clusters chi2
5457 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5458 // in particular, the error associated to the cluster
5459 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5461 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5463 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5464 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5465 // track->SetExtraModule(ilayer,idetExtra);
5467 if(!fPlaneEff->UpDatePlaneEff(found,key))
5468 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5469 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5470 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5471 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5472 Int_t cltype[2]={-999,-999};
5476 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5477 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5480 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5481 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5482 cltype[0]=layer.GetCluster(ci)->GetNy();
5483 cltype[1]=layer.GetCluster(ci)->GetNz();
5485 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5486 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5487 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5488 // It is computed properly by calling the method
5489 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5491 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5492 //if (tmp.PropagateTo(x,0.,0.)) {
5493 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5494 AliCluster c(*layer.GetCluster(ci));
5495 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5496 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5497 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5498 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5499 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5502 fPlaneEff->FillHistos(key,found,tr,clu,cltype);