1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------------------------
18 // Implementation of the ITS tracker class
19 // It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
20 // and fills with them the ESD
21 // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
22 // Current support and development:
23 // Andrea Dainese, andrea.dainese@lnl.infn.it
24 // dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25 // Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26 // Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
27 //-------------------------------------------------------------------------
31 #include <TDatabasePDG.h>
34 #include <TTreeStream.h>
37 #include "AliITSPlaneEff.h"
38 #include "AliITSCalibrationSPD.h"
39 #include "AliITSCalibrationSDD.h"
40 #include "AliITSCalibrationSSD.h"
41 #include "AliCDBEntry.h"
42 #include "AliCDBManager.h"
43 #include "AliAlignObj.h"
44 #include "AliTrackPointArray.h"
45 #include "AliESDVertex.h"
46 #include "AliESDEvent.h"
47 #include "AliESDtrack.h"
50 #include "AliITSChannelStatus.h"
51 #include "AliITSDetTypeRec.h"
52 #include "AliITSRecPoint.h"
53 #include "AliITSgeomTGeo.h"
54 #include "AliITSReconstructor.h"
55 #include "AliITSClusterParam.h"
56 #include "AliITSsegmentation.h"
57 #include "AliITSCalibration.h"
58 #include "AliITSPlaneEffSPD.h"
59 #include "AliITSPlaneEffSDD.h"
60 #include "AliITSPlaneEffSSD.h"
61 #include "AliITSV0Finder.h"
62 #include "AliITStrackerMI.h"
64 ClassImp(AliITStrackerMI)
66 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
68 AliITStrackerMI::AliITStrackerMI():AliTracker(),
78 fLastLayerToTrackTo(0),
81 fTrackingPhase("Default"),
87 fxTimesRhoPipeTrks(0),
88 fxOverX0ShieldTrks(0),
89 fxTimesRhoShieldTrks(0),
91 fxTimesRhoLayerTrks(0),
98 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
99 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
100 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
102 //------------------------------------------------------------------------
103 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
104 fI(AliITSgeomTGeo::GetNLayers()),
113 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
116 fTrackingPhase("Default"),
122 fxTimesRhoPipeTrks(0),
123 fxOverX0ShieldTrks(0),
124 fxTimesRhoShieldTrks(0),
125 fxOverX0LayerTrks(0),
126 fxTimesRhoLayerTrks(0),
128 fITSChannelStatus(0),
131 //--------------------------------------------------------------------
132 //This is the AliITStrackerMI constructor
133 //--------------------------------------------------------------------
135 AliWarning("\"geom\" is actually a dummy argument !");
141 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
142 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
143 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
145 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
146 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
147 Double_t poff=TMath::ATan2(y,x);
149 Double_t r=TMath::Sqrt(x*x + y*y);
151 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
152 r += TMath::Sqrt(x*x + y*y);
153 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
154 r += TMath::Sqrt(x*x + y*y);
155 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
156 r += TMath::Sqrt(x*x + y*y);
159 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
161 for (Int_t j=1; j<nlad+1; j++) {
162 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
163 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
164 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
166 Double_t txyz[3]={0.};
167 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
168 m.LocalToMaster(txyz,xyz);
169 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
170 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
172 if (phi<0) phi+=TMath::TwoPi();
173 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
175 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
176 new(&det) AliITSdetector(r,phi);
177 // compute the real radius (with misalignment)
178 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
180 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
181 mmisal.LocalToMaster(txyz,xyz);
182 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
183 det.SetRmisal(rmisal);
185 } // end loop on detectors
186 } // end loop on ladders
187 fForceSkippingOfLayer[i] = 0;
188 } // end loop on layers
191 fI=AliITSgeomTGeo::GetNLayers();
194 fConstraint[0]=1; fConstraint[1]=0;
196 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
197 AliITSReconstructor::GetRecoParam()->GetYVdef(),
198 AliITSReconstructor::GetRecoParam()->GetZVdef()};
199 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
200 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
201 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
202 SetVertex(xyzVtx,ersVtx);
204 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
205 for (Int_t i=0;i<100000;i++){
206 fBestTrackIndex[i]=0;
209 // store positions of centre of SPD modules (in z)
211 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
212 fSPDdetzcentre[0] = tr[2];
213 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
214 fSPDdetzcentre[1] = tr[2];
215 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
216 fSPDdetzcentre[2] = tr[2];
217 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
218 fSPDdetzcentre[3] = tr[2];
220 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
221 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
222 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
226 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
227 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
229 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
231 // only for plane efficiency evaluation
232 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
233 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
234 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
235 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
236 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
237 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
238 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
239 else fPlaneEff = new AliITSPlaneEffSSD();
240 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
241 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
242 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
245 //------------------------------------------------------------------------
246 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
248 fBestTrack(tracker.fBestTrack),
249 fTrackToFollow(tracker.fTrackToFollow),
250 fTrackHypothesys(tracker.fTrackHypothesys),
251 fBestHypothesys(tracker.fBestHypothesys),
252 fOriginal(tracker.fOriginal),
253 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
254 fPass(tracker.fPass),
255 fAfterV0(tracker.fAfterV0),
256 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
257 fCoefficients(tracker.fCoefficients),
259 fTrackingPhase(tracker.fTrackingPhase),
260 fUseTGeo(tracker.fUseTGeo),
261 fNtracks(tracker.fNtracks),
262 fxOverX0Pipe(tracker.fxOverX0Pipe),
263 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
265 fxTimesRhoPipeTrks(0),
266 fxOverX0ShieldTrks(0),
267 fxTimesRhoShieldTrks(0),
268 fxOverX0LayerTrks(0),
269 fxTimesRhoLayerTrks(0),
270 fDebugStreamer(tracker.fDebugStreamer),
271 fITSChannelStatus(tracker.fITSChannelStatus),
272 fkDetTypeRec(tracker.fkDetTypeRec),
273 fPlaneEff(tracker.fPlaneEff) {
277 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
280 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
281 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
284 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
285 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
288 //------------------------------------------------------------------------
289 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
290 //Assignment operator
291 this->~AliITStrackerMI();
292 new(this) AliITStrackerMI(tracker);
295 //------------------------------------------------------------------------
296 AliITStrackerMI::~AliITStrackerMI()
301 if (fCoefficients) delete [] fCoefficients;
302 DeleteTrksMaterialLUT();
303 if (fDebugStreamer) {
304 //fDebugStreamer->Close();
305 delete fDebugStreamer;
307 if(fITSChannelStatus) delete fITSChannelStatus;
308 if(fPlaneEff) delete fPlaneEff;
310 //------------------------------------------------------------------------
311 void AliITStrackerMI::ReadBadFromDetTypeRec() {
312 //--------------------------------------------------------------------
313 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
315 //--------------------------------------------------------------------
317 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
319 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
321 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
324 if(fITSChannelStatus) delete fITSChannelStatus;
325 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
327 // ITS detectors and chips
328 Int_t i=0,j=0,k=0,ndet=0;
329 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
330 Int_t nBadDetsPerLayer=0;
331 ndet=AliITSgeomTGeo::GetNDetectors(i);
332 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
333 for (k=1; k<ndet+1; k++) {
334 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
335 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
336 if(det.IsBad()) {nBadDetsPerLayer++;}
337 } // end loop on detectors
338 } // end loop on ladders
339 Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
340 } // end loop on layers
344 //------------------------------------------------------------------------
345 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
346 //--------------------------------------------------------------------
347 //This function loads ITS clusters
348 //--------------------------------------------------------------------
349 TBranch *branch=cTree->GetBranch("ITSRecPoints");
351 Error("LoadClusters"," can't get the branch !\n");
355 static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
356 branch->SetAddress(&clusters);
358 Int_t i=0,j=0,ndet=0;
360 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
361 ndet=fgLayers[i].GetNdetectors();
362 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
363 for (; j<jmax; j++) {
364 if (!cTree->GetEvent(j)) continue;
365 Int_t ncl=clusters->GetEntriesFast();
366 SignDeltas(clusters,GetZ());
369 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
370 detector=c->GetDetectorIndex();
372 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
374 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
377 // add dead zone "virtual" cluster in SPD, if there is a cluster within
378 // zwindow cm from the dead zone
379 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
380 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
381 Int_t lab[4] = {0,0,0,detector};
382 Int_t info[3] = {0,0,i};
383 Float_t q = 0.; // this identifies virtual clusters
384 Float_t hit[5] = {xdead,
386 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
387 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
389 Bool_t local = kTRUE;
390 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
391 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
392 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
393 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
394 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
395 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
396 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
397 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
398 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
399 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
400 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
401 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
402 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
403 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
404 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
405 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
406 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
407 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
408 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
410 } // "virtual" clusters in SPD
414 fgLayers[i].ResetRoad(); //road defined by the cluster density
415 fgLayers[i].SortClusters();
418 // check whether we have to skip some layers
419 SetForceSkippingOfLayer();
423 //------------------------------------------------------------------------
424 void AliITStrackerMI::UnloadClusters() {
425 //--------------------------------------------------------------------
426 //This function unloads ITS clusters
427 //--------------------------------------------------------------------
428 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
430 //------------------------------------------------------------------------
431 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
432 //--------------------------------------------------------------------
433 // Publishes all pointers to clusters known to the tracker into the
434 // passed object array.
435 // The ownership is not transfered - the caller is not expected to delete
437 //--------------------------------------------------------------------
439 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
440 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
441 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
448 //------------------------------------------------------------------------
449 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
450 //--------------------------------------------------------------------
451 // Correction for the material between the TPC and the ITS
452 //--------------------------------------------------------------------
453 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
454 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
455 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
456 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
457 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
458 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
459 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
460 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
462 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
468 //------------------------------------------------------------------------
469 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
470 //--------------------------------------------------------------------
471 // This functions reconstructs ITS tracks
472 // The clusters must be already loaded !
473 //--------------------------------------------------------------------
475 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
477 fTrackingPhase="Clusters2Tracks";
479 TObjArray itsTracks(15000);
481 fEsd = event; // store pointer to the esd
483 // temporary (for cosmics)
484 if(event->GetVertex()) {
485 TString title = event->GetVertex()->GetTitle();
486 if(title.Contains("cosmics")) {
487 Double_t xyz[3]={GetX(),GetY(),GetZ()};
488 Double_t exyz[3]={0.1,0.1,0.1};
494 {/* Read ESD tracks */
495 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
496 Int_t nentr=event->GetNumberOfTracks();
497 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
499 AliESDtrack *esd=event->GetTrack(nentr);
500 // ---- for debugging:
501 //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); }
503 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
504 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
505 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
506 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
509 t=new AliITStrackMI(*esd);
510 } catch (const Char_t *msg) {
511 //Warning("Clusters2Tracks",msg);
515 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
516 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
519 // look at the ESD mass hypothesys !
520 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
522 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
524 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
525 //track - can be V0 according to TPC
527 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
531 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
535 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
539 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
544 t->SetReconstructed(kFALSE);
545 itsTracks.AddLast(t);
546 fOriginal.AddLast(t);
548 } /* End Read ESD tracks */
552 Int_t nentr=itsTracks.GetEntriesFast();
553 fTrackHypothesys.Expand(nentr);
554 fBestHypothesys.Expand(nentr);
555 MakeCoefficients(nentr);
556 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
558 // THE TWO TRACKING PASSES
559 for (fPass=0; fPass<2; fPass++) {
560 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
561 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
562 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
563 if (t==0) continue; //this track has been already tracked
564 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
565 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
566 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
567 if (fConstraint[fPass]) {
568 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
569 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
572 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
573 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
575 ResetTrackToFollow(*t);
578 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
581 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
583 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
584 if (!besttrack) continue;
585 besttrack->SetLabel(tpcLabel);
586 // besttrack->CookdEdx();
588 besttrack->SetFakeRatio(1.);
589 CookLabel(besttrack,0.); //For comparison only
590 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
592 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
594 t->SetReconstructed(kTRUE);
596 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
598 GetBestHypothesysMIP(itsTracks);
599 } // end loop on the two tracking passes
601 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
602 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
607 Int_t entries = fTrackHypothesys.GetEntriesFast();
608 for (Int_t ientry=0; ientry<entries; ientry++) {
609 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
610 if (array) array->Delete();
611 delete fTrackHypothesys.RemoveAt(ientry);
614 fTrackHypothesys.Delete();
615 fBestHypothesys.Delete();
617 delete [] fCoefficients;
619 DeleteTrksMaterialLUT();
621 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
623 fTrackingPhase="Default";
627 //------------------------------------------------------------------------
628 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
629 //--------------------------------------------------------------------
630 // This functions propagates reconstructed ITS tracks back
631 // The clusters must be loaded !
632 //--------------------------------------------------------------------
633 fTrackingPhase="PropagateBack";
634 Int_t nentr=event->GetNumberOfTracks();
635 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
638 for (Int_t i=0; i<nentr; i++) {
639 AliESDtrack *esd=event->GetTrack(i);
641 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
642 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
646 t=new AliITStrackMI(*esd);
647 } catch (const Char_t *msg) {
648 //Warning("PropagateBack",msg);
652 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
654 ResetTrackToFollow(*t);
657 // propagate to vertex [SR, GSI 17.02.2003]
658 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
659 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
660 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
661 fTrackToFollow.StartTimeIntegral();
662 // from vertex to outside pipe
663 CorrectForPipeMaterial(&fTrackToFollow,"outward");
665 // Start time integral and add distance from current position to vertex
666 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
667 fTrackToFollow.GetXYZ(xyzTrk);
669 for (Int_t icoord=0; icoord<3; icoord++) {
670 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
673 fTrackToFollow.StartTimeIntegral();
674 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
676 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
677 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
678 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
682 fTrackToFollow.SetLabel(t->GetLabel());
683 //fTrackToFollow.CookdEdx();
684 CookLabel(&fTrackToFollow,0.); //For comparison only
685 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
686 //UseClusters(&fTrackToFollow);
692 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
694 fTrackingPhase="Default";
698 //------------------------------------------------------------------------
699 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
700 //--------------------------------------------------------------------
701 // This functions refits ITS tracks using the
702 // "inward propagated" TPC tracks
703 // The clusters must be loaded !
704 //--------------------------------------------------------------------
705 fTrackingPhase="RefitInward";
707 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
709 Int_t nentr=event->GetNumberOfTracks();
710 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
713 for (Int_t i=0; i<nentr; i++) {
714 AliESDtrack *esd=event->GetTrack(i);
716 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
717 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
718 if (esd->GetStatus()&AliESDtrack::kTPCout)
719 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
723 t=new AliITStrackMI(*esd);
724 } catch (const Char_t *msg) {
725 //Warning("RefitInward",msg);
729 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
730 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
735 ResetTrackToFollow(*t);
736 fTrackToFollow.ResetClusters();
738 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
739 fTrackToFollow.ResetCovariance(10.);
742 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
743 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
745 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
746 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
747 AliDebug(2," refit OK");
748 fTrackToFollow.SetLabel(t->GetLabel());
749 // fTrackToFollow.CookdEdx();
750 CookdEdx(&fTrackToFollow);
752 CookLabel(&fTrackToFollow,0.0); //For comparison only
755 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
756 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
757 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
758 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
759 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
760 Double_t r[3]={0.,0.,0.};
762 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
769 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
771 fTrackingPhase="Default";
775 //------------------------------------------------------------------------
776 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
777 //--------------------------------------------------------------------
778 // Return pointer to a given cluster
779 //--------------------------------------------------------------------
780 Int_t l=(index & 0xf0000000) >> 28;
781 Int_t c=(index & 0x0fffffff) >> 00;
782 return fgLayers[l].GetCluster(c);
784 //------------------------------------------------------------------------
785 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
786 //--------------------------------------------------------------------
787 // Get track space point with index i
788 //--------------------------------------------------------------------
790 Int_t l=(index & 0xf0000000) >> 28;
791 Int_t c=(index & 0x0fffffff) >> 00;
792 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
793 Int_t idet = cl->GetDetectorIndex();
797 cl->GetGlobalXYZ(xyz);
798 cl->GetGlobalCov(cov);
800 p.SetCharge(cl->GetQ());
801 p.SetDriftTime(cl->GetDriftTime());
802 p.SetChargeRatio(cl->GetChargeRatio());
803 p.SetClusterType(cl->GetClusterType());
804 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
807 iLayer = AliGeomManager::kSPD1;
810 iLayer = AliGeomManager::kSPD2;
813 iLayer = AliGeomManager::kSDD1;
816 iLayer = AliGeomManager::kSDD2;
819 iLayer = AliGeomManager::kSSD1;
822 iLayer = AliGeomManager::kSSD2;
825 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
828 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
829 p.SetVolumeID((UShort_t)volid);
832 //------------------------------------------------------------------------
833 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
834 AliTrackPoint& p, const AliESDtrack *t) {
835 //--------------------------------------------------------------------
836 // Get track space point with index i
837 // (assign error estimated during the tracking)
838 //--------------------------------------------------------------------
840 Int_t l=(index & 0xf0000000) >> 28;
841 Int_t c=(index & 0x0fffffff) >> 00;
842 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
843 Int_t idet = cl->GetDetectorIndex();
845 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
847 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
849 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
850 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
851 Double_t alpha = t->GetAlpha();
852 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
853 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
854 phi += alpha-det.GetPhi();
855 Float_t tgphi = TMath::Tan(phi);
857 Float_t tgl = t->GetTgl(); // tgl about const along track
858 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
860 Float_t errtrky,errtrkz,covyz;
861 Bool_t addMisalErr=kFALSE;
862 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
866 cl->GetGlobalXYZ(xyz);
867 // cl->GetGlobalCov(cov);
868 Float_t pos[3] = {0.,0.,0.};
869 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
870 tmpcl.GetGlobalCov(cov);
873 p.SetCharge(cl->GetQ());
874 p.SetDriftTime(cl->GetDriftTime());
875 p.SetChargeRatio(cl->GetChargeRatio());
876 p.SetClusterType(cl->GetClusterType());
878 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
881 iLayer = AliGeomManager::kSPD1;
884 iLayer = AliGeomManager::kSPD2;
887 iLayer = AliGeomManager::kSDD1;
890 iLayer = AliGeomManager::kSDD2;
893 iLayer = AliGeomManager::kSSD1;
896 iLayer = AliGeomManager::kSSD2;
899 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
902 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
904 p.SetVolumeID((UShort_t)volid);
907 //------------------------------------------------------------------------
908 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
910 //--------------------------------------------------------------------
911 // Follow prolongation tree
912 //--------------------------------------------------------------------
914 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
915 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
918 AliESDtrack * esd = otrack->GetESDtrack();
919 if (esd->GetV0Index(0)>0) {
920 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
921 // mapping of ESD track is different as ITS track in Containers
922 // Need something more stable
923 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
924 for (Int_t i=0;i<3;i++){
925 Int_t index = esd->GetV0Index(i);
927 AliESDv0 * vertex = fEsd->GetV0(index);
928 if (vertex->GetStatus()<0) continue; // rejected V0
930 if (esd->GetSign()>0) {
931 vertex->SetIndex(0,esdindex);
933 vertex->SetIndex(1,esdindex);
937 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
939 bestarray = new TObjArray(5);
940 fBestHypothesys.AddAt(bestarray,esdindex);
944 //setup tree of the prolongations
946 static AliITStrackMI tracks[7][100];
947 AliITStrackMI *currenttrack;
948 static AliITStrackMI currenttrack1;
949 static AliITStrackMI currenttrack2;
950 static AliITStrackMI backuptrack;
952 Int_t nindexes[7][100];
953 Float_t normalizedchi2[100];
954 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
955 otrack->SetNSkipped(0);
956 new (&(tracks[6][0])) AliITStrackMI(*otrack);
958 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
959 Int_t modstatus = 1; // found
963 // follow prolongations
964 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
965 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
968 AliITSlayer &layer=fgLayers[ilayer];
969 Double_t r = layer.GetR();
975 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
977 if (ntracks[ilayer]>=100) break;
978 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
979 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
980 if (ntracks[ilayer]>15+ilayer){
981 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
982 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
985 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
987 // material between SSD and SDD, SDD and SPD
989 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
991 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
995 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
996 Int_t idet=layer.FindDetectorIndex(phi,z);
998 Double_t trackGlobXYZ1[3];
999 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1001 // Get the budget to the primary vertex for the current track being prolonged
1002 Double_t budgetToPrimVertex = GetEffectiveThickness();
1004 // check if we allow a prolongation without point
1005 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1007 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1008 // propagate to the layer radius
1009 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1010 if(!vtrack->Propagate(xToGo)) continue;
1011 // apply correction for material of the current layer
1012 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1013 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1014 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1015 vtrack->SetClIndex(ilayer,-1);
1016 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1017 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1018 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1020 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1025 // track outside layer acceptance in z
1026 if (idet<0) continue;
1028 //propagate to the intersection with the detector plane
1029 const AliITSdetector &det=layer.GetDetector(idet);
1030 new(¤ttrack2) AliITStrackMI(currenttrack1);
1031 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1032 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1033 currenttrack1.SetDetectorIndex(idet);
1034 currenttrack2.SetDetectorIndex(idet);
1035 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1038 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1040 // road in global (rphi,z) [i.e. in tracking ref. system]
1041 Double_t zmin,zmax,ymin,ymax;
1042 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1044 // select clusters in road
1045 layer.SelectClusters(zmin,zmax,ymin,ymax);
1046 //********************
1048 // Define criteria for track-cluster association
1049 Double_t msz = currenttrack1.GetSigmaZ2() +
1050 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1051 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1052 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1053 Double_t msy = currenttrack1.GetSigmaY2() +
1054 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1055 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1056 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1058 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1059 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1061 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1062 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1064 msz = 1./msz; // 1/RoadZ^2
1065 msy = 1./msy; // 1/RoadY^2
1069 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1071 const AliITSRecPoint *cl=0;
1073 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1074 Bool_t deadzoneSPD=kFALSE;
1075 currenttrack = ¤ttrack1;
1077 // check if the road contains a dead zone
1078 Bool_t noClusters = kFALSE;
1079 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1080 if (noClusters) AliDebug(2,"no clusters in road");
1081 Double_t dz=0.5*(zmax-zmin);
1082 Double_t dy=0.5*(ymax-ymin);
1083 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1084 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1085 // create a prolongation without clusters (check also if there are no clusters in the road)
1088 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1089 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1090 updatetrack->SetClIndex(ilayer,-1);
1092 modstatus = 5; // no cls in road
1093 } else if (dead==1) {
1094 modstatus = 7; // holes in z in SPD
1095 } else if (dead==2 || dead==3 || dead==4) {
1096 modstatus = 2; // dead from OCDB
1098 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1099 // apply correction for material of the current layer
1100 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1101 if (constrain) { // apply vertex constrain
1102 updatetrack->SetConstrain(constrain);
1103 Bool_t isPrim = kTRUE;
1104 if (ilayer<4) { // check that it's close to the vertex
1105 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1106 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1107 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1108 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1109 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1111 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1113 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1115 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1116 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1118 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1119 updatetrack->SetDeadZoneProbability(ilayer,1.);
1120 } else if (dead==4) { // at least a single dead channel from OCDB
1121 updatetrack->SetDeadZoneProbability(ilayer,0.);
1128 // loop over clusters in the road
1129 while ((cl=layer.GetNextCluster(clidx))!=0) {
1130 if (ntracks[ilayer]>95) break; //space for skipped clusters
1131 Bool_t changedet =kFALSE;
1132 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1133 Int_t idetc=cl->GetDetectorIndex();
1135 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1136 // take into account misalignment (bring track to real detector plane)
1137 Double_t xTrOrig = currenttrack->GetX();
1138 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1139 // a first cut on track-cluster distance
1140 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1141 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1142 { // cluster not associated to track
1143 AliDebug(2,"not associated");
1146 // bring track back to ideal detector plane
1147 if (!currenttrack->Propagate(xTrOrig)) continue;
1148 } else { // have to move track to cluster's detector
1149 const AliITSdetector &detc=layer.GetDetector(idetc);
1150 // a first cut on track-cluster distance
1152 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1153 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1154 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1155 continue; // cluster not associated to track
1157 new (&backuptrack) AliITStrackMI(currenttrack2);
1159 currenttrack =¤ttrack2;
1160 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1161 new (currenttrack) AliITStrackMI(backuptrack);
1165 currenttrack->SetDetectorIndex(idetc);
1166 // Get again the budget to the primary vertex
1167 // for the current track being prolonged, if had to change detector
1168 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1171 // calculate track-clusters chi2
1172 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1174 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1175 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1176 if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1177 if (ntracks[ilayer]>=100) continue;
1178 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1179 updatetrack->SetClIndex(ilayer,-1);
1180 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1182 if (cl->GetQ()!=0) { // real cluster
1183 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1184 AliDebug(2,"update failed");
1187 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1188 modstatus = 1; // found
1189 } else { // virtual cluster in dead zone
1190 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1191 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1192 modstatus = 7; // holes in z in SPD
1196 Float_t xlocnewdet,zlocnewdet;
1197 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1198 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1201 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1203 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1205 // apply correction for material of the current layer
1206 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1208 if (constrain) { // apply vertex constrain
1209 updatetrack->SetConstrain(constrain);
1210 Bool_t isPrim = kTRUE;
1211 if (ilayer<4) { // check that it's close to the vertex
1212 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1213 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1214 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1215 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1216 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1218 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1219 } //apply vertex constrain
1221 } // create new hypothesis
1223 AliDebug(2,"chi2 too large");
1226 } // loop over possible prolongations
1228 // allow one prolongation without clusters
1229 if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1230 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1231 // apply correction for material of the current layer
1232 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1233 vtrack->SetClIndex(ilayer,-1);
1234 modstatus = 3; // skipped
1235 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1236 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1237 vtrack->IncrementNSkipped();
1241 // allow one prolongation without clusters for tracks with |tgl|>1.1
1242 if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) { //big theta - for low flux
1243 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1244 // apply correction for material of the current layer
1245 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1246 vtrack->SetClIndex(ilayer,-1);
1247 modstatus = 3; // skipped
1248 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1249 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1250 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1255 } // loop over tracks in layer ilayer+1
1257 //loop over track candidates for the current layer
1263 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1264 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1265 if (normalizedchi2[itrack] <
1266 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1270 if (constrain) { // constrain
1271 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1273 } else { // !constrain
1274 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1279 // sort tracks by increasing normalized chi2
1280 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1281 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1282 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1283 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1284 } // end loop over layers
1288 // Now select tracks to be kept
1290 Int_t max = constrain ? 20 : 5;
1292 // tracks that reach layer 0 (SPD inner)
1293 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1294 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1295 if (track.GetNumberOfClusters()<2) continue;
1296 if (!constrain && track.GetNormChi2(0) >
1297 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1300 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1303 // tracks that reach layer 1 (SPD outer)
1304 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1305 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1306 if (track.GetNumberOfClusters()<4) continue;
1307 if (!constrain && track.GetNormChi2(1) >
1308 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1309 if (constrain) track.IncrementNSkipped();
1311 track.SetD(0,track.GetD(GetX(),GetY()));
1312 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1313 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1314 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1317 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1320 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1322 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1323 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1324 if (track.GetNumberOfClusters()<3) continue;
1325 if (!constrain && track.GetNormChi2(2) >
1326 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1327 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1329 track.SetD(0,track.GetD(GetX(),GetY()));
1330 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1331 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1332 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1335 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1341 // register best track of each layer - important for V0 finder
1343 for (Int_t ilayer=0;ilayer<5;ilayer++){
1344 if (ntracks[ilayer]==0) continue;
1345 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1346 if (track.GetNumberOfClusters()<1) continue;
1347 CookLabel(&track,0);
1348 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1352 // update TPC V0 information
1354 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1355 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1356 for (Int_t i=0;i<3;i++){
1357 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1358 if (index==0) break;
1359 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1360 if (vertex->GetStatus()<0) continue; // rejected V0
1362 if (otrack->GetSign()>0) {
1363 vertex->SetIndex(0,esdindex);
1366 vertex->SetIndex(1,esdindex);
1368 //find nearest layer with track info
1369 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1370 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1371 Int_t nearest = nearestold;
1372 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1373 if (ntracks[nearest]==0){
1378 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1379 if (nearestold<5&&nearest<5){
1380 Bool_t accept = track.GetNormChi2(nearest)<10;
1382 if (track.GetSign()>0) {
1383 vertex->SetParamP(track);
1384 vertex->Update(fprimvertex);
1385 //vertex->SetIndex(0,track.fESDtrack->GetID());
1386 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1388 vertex->SetParamN(track);
1389 vertex->Update(fprimvertex);
1390 //vertex->SetIndex(1,track.fESDtrack->GetID());
1391 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1393 vertex->SetStatus(vertex->GetStatus()+1);
1395 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1402 //------------------------------------------------------------------------
1403 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1405 //--------------------------------------------------------------------
1408 return fgLayers[layer];
1410 //------------------------------------------------------------------------
1411 AliITStrackerMI::AliITSlayer::AliITSlayer():
1441 //--------------------------------------------------------------------
1442 //default AliITSlayer constructor
1443 //--------------------------------------------------------------------
1444 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1445 fClusterWeight[i]=0;
1446 fClusterTracks[0][i]=-1;
1447 fClusterTracks[1][i]=-1;
1448 fClusterTracks[2][i]=-1;
1449 fClusterTracks[3][i]=-1;
1452 //------------------------------------------------------------------------
1453 AliITStrackerMI::AliITSlayer::
1454 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1483 //--------------------------------------------------------------------
1484 //main AliITSlayer constructor
1485 //--------------------------------------------------------------------
1486 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1487 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1489 //------------------------------------------------------------------------
1490 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1492 fPhiOffset(layer.fPhiOffset),
1493 fNladders(layer.fNladders),
1494 fZOffset(layer.fZOffset),
1495 fNdetectors(layer.fNdetectors),
1496 fDetectors(layer.fDetectors),
1501 fClustersCs(layer.fClustersCs),
1502 fClusterIndexCs(layer.fClusterIndexCs),
1506 fCurrentSlice(layer.fCurrentSlice),
1514 fAccepted(layer.fAccepted),
1516 fMaxSigmaClY(layer.fMaxSigmaClY),
1517 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1518 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1522 //------------------------------------------------------------------------
1523 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1524 //--------------------------------------------------------------------
1525 // AliITSlayer destructor
1526 //--------------------------------------------------------------------
1527 delete [] fDetectors;
1528 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1529 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1530 fClusterWeight[i]=0;
1531 fClusterTracks[0][i]=-1;
1532 fClusterTracks[1][i]=-1;
1533 fClusterTracks[2][i]=-1;
1534 fClusterTracks[3][i]=-1;
1537 //------------------------------------------------------------------------
1538 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1539 //--------------------------------------------------------------------
1540 // This function removes loaded clusters
1541 //--------------------------------------------------------------------
1542 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1543 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1544 fClusterWeight[i]=0;
1545 fClusterTracks[0][i]=-1;
1546 fClusterTracks[1][i]=-1;
1547 fClusterTracks[2][i]=-1;
1548 fClusterTracks[3][i]=-1;
1554 //------------------------------------------------------------------------
1555 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1556 //--------------------------------------------------------------------
1557 // This function reset weights of the clusters
1558 //--------------------------------------------------------------------
1559 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1560 fClusterWeight[i]=0;
1561 fClusterTracks[0][i]=-1;
1562 fClusterTracks[1][i]=-1;
1563 fClusterTracks[2][i]=-1;
1564 fClusterTracks[3][i]=-1;
1566 for (Int_t i=0; i<fN;i++) {
1567 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1568 if (cl&&cl->IsUsed()) cl->Use();
1572 //------------------------------------------------------------------------
1573 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1574 //--------------------------------------------------------------------
1575 // This function calculates the road defined by the cluster density
1576 //--------------------------------------------------------------------
1578 for (Int_t i=0; i<fN; i++) {
1579 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1581 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1583 //------------------------------------------------------------------------
1584 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1585 //--------------------------------------------------------------------
1586 //This function adds a cluster to this layer
1587 //--------------------------------------------------------------------
1588 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1589 ::Error("InsertCluster","Too many clusters !\n");
1595 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1597 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1598 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1599 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1600 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1601 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1602 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1605 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1606 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1607 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1608 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1612 //------------------------------------------------------------------------
1613 void AliITStrackerMI::AliITSlayer::SortClusters()
1618 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1619 Float_t *z = new Float_t[fN];
1620 Int_t * index = new Int_t[fN];
1622 fMaxSigmaClY=0.; //AD
1623 fMaxSigmaClZ=0.; //AD
1625 for (Int_t i=0;i<fN;i++){
1626 z[i] = fClusters[i]->GetZ();
1627 // save largest errors in y and z for this layer
1628 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1629 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1631 TMath::Sort(fN,z,index,kFALSE);
1632 for (Int_t i=0;i<fN;i++){
1633 clusters[i] = fClusters[index[i]];
1636 for (Int_t i=0;i<fN;i++){
1637 fClusters[i] = clusters[i];
1638 fZ[i] = fClusters[i]->GetZ();
1639 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1640 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1641 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1651 for (Int_t i=0;i<fN;i++){
1652 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1653 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1654 fClusterIndex[i] = i;
1658 fDy5 = (fYB[1]-fYB[0])/5.;
1659 fDy10 = (fYB[1]-fYB[0])/10.;
1660 fDy20 = (fYB[1]-fYB[0])/20.;
1661 for (Int_t i=0;i<6;i++) fN5[i] =0;
1662 for (Int_t i=0;i<11;i++) fN10[i]=0;
1663 for (Int_t i=0;i<21;i++) fN20[i]=0;
1665 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;}
1666 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;}
1667 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;}
1670 for (Int_t i=0;i<fN;i++)
1671 for (Int_t irot=-1;irot<=1;irot++){
1672 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1674 for (Int_t slice=0; slice<6;slice++){
1675 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1676 fClusters5[slice][fN5[slice]] = fClusters[i];
1677 fY5[slice][fN5[slice]] = curY;
1678 fZ5[slice][fN5[slice]] = fZ[i];
1679 fClusterIndex5[slice][fN5[slice]]=i;
1684 for (Int_t slice=0; slice<11;slice++){
1685 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1686 fClusters10[slice][fN10[slice]] = fClusters[i];
1687 fY10[slice][fN10[slice]] = curY;
1688 fZ10[slice][fN10[slice]] = fZ[i];
1689 fClusterIndex10[slice][fN10[slice]]=i;
1694 for (Int_t slice=0; slice<21;slice++){
1695 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1696 fClusters20[slice][fN20[slice]] = fClusters[i];
1697 fY20[slice][fN20[slice]] = curY;
1698 fZ20[slice][fN20[slice]] = fZ[i];
1699 fClusterIndex20[slice][fN20[slice]]=i;
1706 // consistency check
1708 for (Int_t i=0;i<fN-1;i++){
1714 for (Int_t slice=0;slice<21;slice++)
1715 for (Int_t i=0;i<fN20[slice]-1;i++){
1716 if (fZ20[slice][i]>fZ20[slice][i+1]){
1723 //------------------------------------------------------------------------
1724 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1725 //--------------------------------------------------------------------
1726 // This function returns the index of the nearest cluster
1727 //--------------------------------------------------------------------
1730 if (fCurrentSlice<0) {
1739 if (ncl==0) return 0;
1740 Int_t b=0, e=ncl-1, m=(b+e)/2;
1741 for (; b<e; m=(b+e)/2) {
1742 // if (z > fClusters[m]->GetZ()) b=m+1;
1743 if (z > zcl[m]) b=m+1;
1748 //------------------------------------------------------------------------
1749 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 {
1750 //--------------------------------------------------------------------
1751 // This function computes the rectangular road for this track
1752 //--------------------------------------------------------------------
1755 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1756 // take into account the misalignment: propagate track to misaligned detector plane
1757 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1759 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1760 TMath::Sqrt(track->GetSigmaZ2() +
1761 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1762 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1763 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1764 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1765 TMath::Sqrt(track->GetSigmaY2() +
1766 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1767 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1768 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1770 // track at boundary between detectors, enlarge road
1771 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1772 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1773 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1774 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1775 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1776 Float_t tgl = TMath::Abs(track->GetTgl());
1777 if (tgl > 1.) tgl=1.;
1778 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1779 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1780 Float_t snp = TMath::Abs(track->GetSnp());
1781 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1782 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1785 // add to the road a term (up to 2-3 mm) to deal with misalignments
1786 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1787 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1789 Double_t r = fgLayers[ilayer].GetR();
1790 zmin = track->GetZ() - dz;
1791 zmax = track->GetZ() + dz;
1792 ymin = track->GetY() + r*det.GetPhi() - dy;
1793 ymax = track->GetY() + r*det.GetPhi() + dy;
1795 // bring track back to idead detector plane
1796 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1800 //------------------------------------------------------------------------
1801 void AliITStrackerMI::AliITSlayer::
1802 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1803 //--------------------------------------------------------------------
1804 // This function sets the "window"
1805 //--------------------------------------------------------------------
1807 Double_t circle=2*TMath::Pi()*fR;
1813 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1814 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1815 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1816 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1817 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1819 Float_t ymiddle = (fYmin+fYmax)*0.5;
1820 if (ymiddle<fYB[0]) {
1821 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1822 } else if (ymiddle>fYB[1]) {
1823 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1829 fClustersCs = fClusters;
1830 fClusterIndexCs = fClusterIndex;
1836 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1837 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1838 if (slice<0) slice=0;
1839 if (slice>20) slice=20;
1840 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1842 fCurrentSlice=slice;
1843 fClustersCs = fClusters20[fCurrentSlice];
1844 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1845 fYcs = fY20[fCurrentSlice];
1846 fZcs = fZ20[fCurrentSlice];
1847 fNcs = fN20[fCurrentSlice];
1852 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1853 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1854 if (slice<0) slice=0;
1855 if (slice>10) slice=10;
1856 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1858 fCurrentSlice=slice;
1859 fClustersCs = fClusters10[fCurrentSlice];
1860 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1861 fYcs = fY10[fCurrentSlice];
1862 fZcs = fZ10[fCurrentSlice];
1863 fNcs = fN10[fCurrentSlice];
1868 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1869 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1870 if (slice<0) slice=0;
1871 if (slice>5) slice=5;
1872 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1874 fCurrentSlice=slice;
1875 fClustersCs = fClusters5[fCurrentSlice];
1876 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1877 fYcs = fY5[fCurrentSlice];
1878 fZcs = fZ5[fCurrentSlice];
1879 fNcs = fN5[fCurrentSlice];
1883 fI = FindClusterIndex(fZmin);
1884 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
1890 //------------------------------------------------------------------------
1891 Int_t AliITStrackerMI::AliITSlayer::
1892 FindDetectorIndex(Double_t phi, Double_t z) const {
1893 //--------------------------------------------------------------------
1894 //This function finds the detector crossed by the track
1895 //--------------------------------------------------------------------
1897 if (fZOffset<0) // old geometry
1898 dphi = -(phi-fPhiOffset);
1899 else // new geometry
1900 dphi = phi-fPhiOffset;
1903 if (dphi < 0) dphi += 2*TMath::Pi();
1904 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1905 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1906 if (np>=fNladders) np-=fNladders;
1907 if (np<0) np+=fNladders;
1910 Double_t dz=fZOffset-z;
1911 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1912 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1913 if (nz>=fNdetectors) return -1;
1914 if (nz<0) return -1;
1916 // ad hoc correction for 3rd ladder of SDD inner layer,
1917 // which is reversed (rotated by pi around local y)
1918 // this correction is OK only from AliITSv11Hybrid onwards
1919 if (GetR()>12. && GetR()<20.) { // SDD inner
1920 if(np==2) { // 3rd ladder
1921 nz = (fNdetectors-1) - nz;
1924 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1927 return np*fNdetectors + nz;
1929 //------------------------------------------------------------------------
1930 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1932 //--------------------------------------------------------------------
1933 // This function returns clusters within the "window"
1934 //--------------------------------------------------------------------
1936 if (fCurrentSlice<0) {
1937 Double_t rpi2 = 2.*fR*TMath::Pi();
1938 for (Int_t i=fI; i<fImax; i++) {
1941 if (fYmax<y) y -= rpi2;
1942 if (fYmin>y) y += rpi2;
1943 if (y<fYmin) continue;
1944 if (y>fYmax) continue;
1946 // skip clusters that are in "extended" road but they
1947 // 3sigma error does not touch the original road
1948 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
1949 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
1951 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1954 return fClusters[i];
1957 for (Int_t i=fI; i<fImax; i++) {
1958 if (fYcs[i]<fYmin) continue;
1959 if (fYcs[i]>fYmax) continue;
1960 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1961 ci=fClusterIndexCs[i];
1963 return fClustersCs[i];
1968 //------------------------------------------------------------------------
1969 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1971 //--------------------------------------------------------------------
1972 // This function returns the layer thickness at this point (units X0)
1973 //--------------------------------------------------------------------
1975 x0=AliITSRecoParam::GetX0Air();
1976 if (43<fR&&fR<45) { //SSD2
1979 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1980 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1981 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1982 for (Int_t i=0; i<12; i++) {
1983 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1984 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1988 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1989 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1993 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1994 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1997 if (37<fR&&fR<41) { //SSD1
2000 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2001 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2002 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2003 for (Int_t i=0; i<11; i++) {
2004 if (TMath::Abs(z-3.9*i)<0.15) {
2005 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2009 if (TMath::Abs(z+3.9*i)<0.15) {
2010 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2014 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2015 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2018 if (13<fR&&fR<26) { //SDD
2021 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2023 if (TMath::Abs(y-1.80)<0.55) {
2025 for (Int_t j=0; j<20; j++) {
2026 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2027 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2030 if (TMath::Abs(y+1.80)<0.55) {
2032 for (Int_t j=0; j<20; j++) {
2033 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2034 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2038 for (Int_t i=0; i<4; i++) {
2039 if (TMath::Abs(z-7.3*i)<0.60) {
2041 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2044 if (TMath::Abs(z+7.3*i)<0.60) {
2046 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2051 if (6<fR&&fR<8) { //SPD2
2052 Double_t dd=0.0063; x0=21.5;
2054 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2055 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2057 if (3<fR&&fR<5) { //SPD1
2058 Double_t dd=0.0063; x0=21.5;
2060 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2061 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2066 //------------------------------------------------------------------------
2067 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2069 fRmisal(det.fRmisal),
2071 fSinPhi(det.fSinPhi),
2072 fCosPhi(det.fCosPhi),
2078 fNChips(det.fNChips),
2079 fChipIsBad(det.fChipIsBad)
2083 //------------------------------------------------------------------------
2084 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2085 const AliITSDetTypeRec *detTypeRec)
2087 //--------------------------------------------------------------------
2088 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2089 //--------------------------------------------------------------------
2091 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2092 // while in the tracker they start from 0 for each layer
2093 for(Int_t il=0; il<ilayer; il++)
2094 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2097 if (ilayer==0 || ilayer==1) { // ---------- SPD
2099 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2101 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2104 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2108 // Get calibration from AliITSDetTypeRec
2109 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2110 calib->SetModuleIndex(idet);
2111 AliITSCalibration *calibSPDdead = 0;
2112 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2113 if (calib->IsBad() ||
2114 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2117 // printf("lay %d bad %d\n",ilayer,idet);
2120 // Get segmentation from AliITSDetTypeRec
2121 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2123 // Read info about bad chips
2124 fNChips = segm->GetMaximumChipIndex()+1;
2125 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2126 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2127 fChipIsBad = new Bool_t[fNChips];
2128 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2129 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2130 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2131 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2136 //------------------------------------------------------------------------
2137 Double_t AliITStrackerMI::GetEffectiveThickness()
2139 //--------------------------------------------------------------------
2140 // Returns the thickness between the current layer and the vertex (units X0)
2141 //--------------------------------------------------------------------
2144 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2145 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2146 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2150 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2151 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2155 Double_t xn=fgLayers[fI].GetR();
2156 for (Int_t i=0; i<fI; i++) {
2157 Double_t xi=fgLayers[i].GetR();
2158 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2164 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2165 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2168 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2169 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2173 //------------------------------------------------------------------------
2174 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2175 //-------------------------------------------------------------------
2176 // This function returns number of clusters within the "window"
2177 //--------------------------------------------------------------------
2179 for (Int_t i=fI; i<fN; i++) {
2180 const AliITSRecPoint *c=fClusters[i];
2181 if (c->GetZ() > fZmax) break;
2182 if (c->IsUsed()) continue;
2183 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2184 Double_t y=fR*det.GetPhi() + c->GetY();
2186 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2187 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2189 if (y<fYmin) continue;
2190 if (y>fYmax) continue;
2195 //------------------------------------------------------------------------
2196 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2197 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2199 //--------------------------------------------------------------------
2200 // This function refits the track "track" at the position "x" using
2201 // the clusters from "clusters"
2202 // If "extra"==kTRUE,
2203 // the clusters from overlapped modules get attached to "track"
2204 // If "planeff"==kTRUE,
2205 // special approach for plane efficiency evaluation is applyed
2206 //--------------------------------------------------------------------
2208 Int_t index[AliITSgeomTGeo::kNLayers];
2210 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2211 Int_t nc=clusters->GetNumberOfClusters();
2212 for (k=0; k<nc; k++) {
2213 Int_t idx=clusters->GetClusterIndex(k);
2214 Int_t ilayer=(idx&0xf0000000)>>28;
2218 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2220 //------------------------------------------------------------------------
2221 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2222 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2224 //--------------------------------------------------------------------
2225 // This function refits the track "track" at the position "x" using
2226 // the clusters from array
2227 // If "extra"==kTRUE,
2228 // the clusters from overlapped modules get attached to "track"
2229 // If "planeff"==kTRUE,
2230 // special approach for plane efficiency evaluation is applyed
2231 //--------------------------------------------------------------------
2232 Int_t index[AliITSgeomTGeo::kNLayers];
2234 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2236 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2237 index[k]=clusters[k];
2240 // special for cosmics: check which the innermost layer crossed
2242 Int_t innermostlayer=5;
2243 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2244 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2245 if(drphi < fgLayers[innermostlayer].GetR()) break;
2247 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2249 Int_t modstatus=1; // found
2251 Int_t from, to, step;
2252 if (xx > track->GetX()) {
2253 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2256 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2259 TString dir = (step>0 ? "outward" : "inward");
2261 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2262 AliITSlayer &layer=fgLayers[ilayer];
2263 Double_t r=layer.GetR();
2265 if (step<0 && xx>r) break;
2267 // material between SSD and SDD, SDD and SPD
2268 Double_t hI=ilayer-0.5*step;
2269 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2270 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2271 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2272 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2275 Double_t oldGlobXYZ[3];
2276 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2278 // continue if we are already beyond this layer
2279 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2280 if(step>0 && oldGlobR > r) continue; // going outward
2281 if(step<0 && oldGlobR < r) continue; // going inward
2284 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2286 Int_t idet=layer.FindDetectorIndex(phi,z);
2288 // check if we allow a prolongation without point for large-eta tracks
2289 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2291 modstatus = 4; // out in z
2292 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2293 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2296 // apply correction for material of the current layer
2297 // add time if going outward
2298 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2302 if (idet<0) return kFALSE;
2305 const AliITSdetector &det=layer.GetDetector(idet);
2306 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2308 track->SetDetectorIndex(idet);
2309 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2311 Double_t dz,zmin,zmax,dy,ymin,ymax;
2313 const AliITSRecPoint *clAcc=0;
2314 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2316 Int_t idx=index[ilayer];
2317 if (idx>=0) { // cluster in this layer
2318 modstatus = 6; // no refit
2319 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2321 if (idet != cl->GetDetectorIndex()) {
2322 idet=cl->GetDetectorIndex();
2323 const AliITSdetector &detc=layer.GetDetector(idet);
2324 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2325 track->SetDetectorIndex(idet);
2326 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2328 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2329 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2333 modstatus = 1; // found
2338 } else { // no cluster in this layer
2340 modstatus = 3; // skipped
2341 // Plane Eff determination:
2342 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2343 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2344 UseTrackForPlaneEff(track,ilayer);
2347 modstatus = 5; // no cls in road
2349 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2350 dz = 0.5*(zmax-zmin);
2351 dy = 0.5*(ymax-ymin);
2352 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2353 if (dead==1) modstatus = 7; // holes in z in SPD
2354 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2359 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2360 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2362 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2365 if (extra) { // search for extra clusters in overlapped modules
2366 AliITStrackV2 tmp(*track);
2367 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2368 layer.SelectClusters(zmin,zmax,ymin,ymax);
2370 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2372 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2373 Double_t tolerance=0.1;
2374 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2375 // only clusters in another module! (overlaps)
2376 idetExtra = clExtra->GetDetectorIndex();
2377 if (idet == idetExtra) continue;
2379 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2381 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2382 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2383 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2384 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2386 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2387 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2390 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2391 track->SetExtraModule(ilayer,idetExtra);
2393 } // end search for extra clusters in overlapped modules
2395 // Correct for material of the current layer
2397 // add time if going outward
2398 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2400 } // end loop on layers
2402 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2406 //------------------------------------------------------------------------
2407 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2410 // calculate normalized chi2
2411 // return NormalizedChi2(track,0);
2414 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2415 // track->fdEdxMismatch=0;
2416 Float_t dedxmismatch =0;
2417 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2419 for (Int_t i = 0;i<6;i++){
2420 if (track->GetClIndex(i)>=0){
2421 Float_t cerry, cerrz;
2422 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2424 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2427 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2428 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2429 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2431 cchi2+=(0.5-ratio)*10.;
2432 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2433 dedxmismatch+=(0.5-ratio)*10.;
2437 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2438 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2439 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2440 if (i<2) chi2+=2*cl->GetDeltaProbability();
2446 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2447 track->SetdEdxMismatch(dedxmismatch);
2451 for (Int_t i = 0;i<4;i++){
2452 if (track->GetClIndex(i)>=0){
2453 Float_t cerry, cerrz;
2454 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2455 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2458 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2459 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2463 for (Int_t i = 4;i<6;i++){
2464 if (track->GetClIndex(i)>=0){
2465 Float_t cerry, cerrz;
2466 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2467 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2470 Float_t cerryb, cerrzb;
2471 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2472 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2475 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2476 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2481 if (track->GetESDtrack()->GetTPCsignal()>85){
2482 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2484 chi2+=(0.5-ratio)*5.;
2487 chi2+=(ratio-2.0)*3;
2491 Double_t match = TMath::Sqrt(track->GetChi22());
2492 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2493 if (!track->GetConstrain()) {
2494 if (track->GetNumberOfClusters()>2) {
2495 match/=track->GetNumberOfClusters()-2.;
2500 if (match<0) match=0;
2502 // penalty factor for missing points (NDeadZone>0), but no penalty
2503 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2504 // or there is a dead from OCDB)
2505 Float_t deadzonefactor = 0.;
2506 if (track->GetNDeadZone()>0.) {
2507 Float_t sumDeadZoneProbability=0;
2508 for(Int_t ilay=0;ilay<6;ilay++) sumDeadZoneProbability+=track->GetDeadZoneProbability(ilay);
2509 Float_t nDeadZoneWithProbNot1=(Float_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2510 if(nDeadZoneWithProbNot1>0.) {
2511 Float_t deadZoneProbability = sumDeadZoneProbability - (Float_t)((Int_t)sumDeadZoneProbability);
2512 deadZoneProbability /= nDeadZoneWithProbNot1;
2513 deadzonefactor = 3.*(1.1-deadZoneProbability);
2517 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2518 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2519 1./(1.+track->GetNSkipped()));
2523 //------------------------------------------------------------------------
2524 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2527 // return matching chi2 between two tracks
2528 Double_t largeChi2=1000.;
2530 AliITStrackMI track3(*track2);
2531 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2533 vec(0,0)=track1->GetY() - track3.GetY();
2534 vec(1,0)=track1->GetZ() - track3.GetZ();
2535 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2536 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2537 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2540 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2541 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2542 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2543 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2544 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2546 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2547 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2548 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2549 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2551 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2552 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2553 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2555 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2556 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2558 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2561 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2562 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2565 //------------------------------------------------------------------------
2566 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2569 // return probability that given point (characterized by z position and error)
2570 // is in SPD dead zone
2572 Double_t probability = 0.;
2573 Double_t absz = TMath::Abs(zpos);
2574 Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2575 0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2576 if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2577 Double_t zmin, zmax;
2578 if (zpos<-6.) { // dead zone at z = -7
2579 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2580 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2581 } else if (zpos>6.) { // dead zone at z = +7
2582 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2583 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2584 } else if (absz<2.) { // dead zone at z = 0
2585 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2586 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2591 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2593 probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2594 TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2597 //------------------------------------------------------------------------
2598 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2601 // calculate normalized chi2
2603 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2605 for (Int_t i = 0;i<6;i++){
2606 if (TMath::Abs(track->GetDy(i))>0){
2607 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2608 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2611 else{chi2[i]=10000;}
2614 TMath::Sort(6,chi2,index,kFALSE);
2615 Float_t max = float(ncl)*fac-1.;
2616 Float_t sumchi=0, sumweight=0;
2617 for (Int_t i=0;i<max+1;i++){
2618 Float_t weight = (i<max)?1.:(max+1.-i);
2619 sumchi+=weight*chi2[index[i]];
2622 Double_t normchi2 = sumchi/sumweight;
2625 //------------------------------------------------------------------------
2626 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2629 // calculate normalized chi2
2630 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2633 for (Int_t i=0;i<6;i++){
2634 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2635 Double_t sy1 = forwardtrack->GetSigmaY(i);
2636 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2637 Double_t sy2 = backtrack->GetSigmaY(i);
2638 Double_t sz2 = backtrack->GetSigmaZ(i);
2639 if (i<2){ sy2=1000.;sz2=1000;}
2641 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2642 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2644 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2645 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2647 res+= nz0*nz0+ny0*ny0;
2650 if (npoints>1) return
2651 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2652 //2*forwardtrack->fNUsed+
2653 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2654 1./(1.+forwardtrack->GetNSkipped()));
2657 //------------------------------------------------------------------------
2658 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2659 //--------------------------------------------------------------------
2660 // Return pointer to a given cluster
2661 //--------------------------------------------------------------------
2662 Int_t l=(index & 0xf0000000) >> 28;
2663 Int_t c=(index & 0x0fffffff) >> 00;
2664 return fgLayers[l].GetWeight(c);
2666 //------------------------------------------------------------------------
2667 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2669 //---------------------------------------------
2670 // register track to the list
2672 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2675 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2676 Int_t index = track->GetClusterIndex(icluster);
2677 Int_t l=(index & 0xf0000000) >> 28;
2678 Int_t c=(index & 0x0fffffff) >> 00;
2679 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2680 for (Int_t itrack=0;itrack<4;itrack++){
2681 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2682 fgLayers[l].SetClusterTracks(itrack,c,id);
2688 //------------------------------------------------------------------------
2689 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2691 //---------------------------------------------
2692 // unregister track from the list
2693 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2694 Int_t index = track->GetClusterIndex(icluster);
2695 Int_t l=(index & 0xf0000000) >> 28;
2696 Int_t c=(index & 0x0fffffff) >> 00;
2697 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2698 for (Int_t itrack=0;itrack<4;itrack++){
2699 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2700 fgLayers[l].SetClusterTracks(itrack,c,-1);
2705 //------------------------------------------------------------------------
2706 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2708 //-------------------------------------------------------------
2709 //get number of shared clusters
2710 //-------------------------------------------------------------
2712 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2713 // mean number of clusters
2714 Float_t *ny = GetNy(id), *nz = GetNz(id);
2717 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2718 Int_t index = track->GetClusterIndex(icluster);
2719 Int_t l=(index & 0xf0000000) >> 28;
2720 Int_t c=(index & 0x0fffffff) >> 00;
2721 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2723 printf("problem\n");
2725 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2729 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2730 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2731 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2733 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2736 deltan = (cl->GetNz()-nz[l]);
2738 if (deltan>2.0) continue; // extended - highly probable shared cluster
2739 weight = 2./TMath::Max(3.+deltan,2.);
2741 for (Int_t itrack=0;itrack<4;itrack++){
2742 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2744 clist[l] = (AliITSRecPoint*)GetCluster(index);
2750 track->SetNUsed(shared);
2753 //------------------------------------------------------------------------
2754 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2757 // find first shared track
2759 // mean number of clusters
2760 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2762 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2763 Int_t sharedtrack=100000;
2764 Int_t tracks[24],trackindex=0;
2765 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2767 for (Int_t icluster=0;icluster<6;icluster++){
2768 if (clusterlist[icluster]<0) continue;
2769 Int_t index = clusterlist[icluster];
2770 Int_t l=(index & 0xf0000000) >> 28;
2771 Int_t c=(index & 0x0fffffff) >> 00;
2773 printf("problem\n");
2775 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2776 //if (l>3) continue;
2777 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2780 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2781 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2782 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2784 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2787 deltan = (cl->GetNz()-nz[l]);
2789 if (deltan>2.0) continue; // extended - highly probable shared cluster
2791 for (Int_t itrack=3;itrack>=0;itrack--){
2792 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2793 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2794 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2799 if (trackindex==0) return -1;
2801 sharedtrack = tracks[0];
2803 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2806 Int_t tracks2[24], cluster[24];
2807 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2810 for (Int_t i=0;i<trackindex;i++){
2811 if (tracks[i]<0) continue;
2812 tracks2[index] = tracks[i];
2814 for (Int_t j=i+1;j<trackindex;j++){
2815 if (tracks[j]<0) continue;
2816 if (tracks[j]==tracks[i]){
2824 for (Int_t i=0;i<index;i++){
2825 if (cluster[index]>max) {
2826 sharedtrack=tracks2[index];
2833 if (sharedtrack>=100000) return -1;
2835 // make list of overlaps
2837 for (Int_t icluster=0;icluster<6;icluster++){
2838 if (clusterlist[icluster]<0) continue;
2839 Int_t index = clusterlist[icluster];
2840 Int_t l=(index & 0xf0000000) >> 28;
2841 Int_t c=(index & 0x0fffffff) >> 00;
2842 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2843 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2845 if (cl->GetNy()>2) continue;
2846 if (cl->GetNz()>2) continue;
2849 if (cl->GetNy()>3) continue;
2850 if (cl->GetNz()>3) continue;
2853 for (Int_t itrack=3;itrack>=0;itrack--){
2854 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2855 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2863 //------------------------------------------------------------------------
2864 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2866 // try to find track hypothesys without conflicts
2867 // with minimal chi2;
2868 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2869 Int_t entries1 = arr1->GetEntriesFast();
2870 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2871 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2872 Int_t entries2 = arr2->GetEntriesFast();
2873 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2875 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2876 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2877 if (track10->Pt()>0.5+track20->Pt()) return track10;
2879 for (Int_t itrack=0;itrack<entries1;itrack++){
2880 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2881 UnRegisterClusterTracks(track,trackID1);
2884 for (Int_t itrack=0;itrack<entries2;itrack++){
2885 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2886 UnRegisterClusterTracks(track,trackID2);
2890 Float_t maxconflicts=6;
2891 Double_t maxchi2 =1000.;
2893 // get weight of hypothesys - using best hypothesy
2896 Int_t list1[6],list2[6];
2897 AliITSRecPoint *clist1[6], *clist2[6] ;
2898 RegisterClusterTracks(track10,trackID1);
2899 RegisterClusterTracks(track20,trackID2);
2900 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2901 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2902 UnRegisterClusterTracks(track10,trackID1);
2903 UnRegisterClusterTracks(track20,trackID2);
2906 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2907 Float_t nerry[6],nerrz[6];
2908 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2909 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2910 for (Int_t i=0;i<6;i++){
2911 if ( (erry1[i]>0) && (erry2[i]>0)) {
2912 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2913 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2915 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2916 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2918 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2919 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2920 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2923 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2924 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2925 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2933 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2934 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2935 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2936 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2938 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2939 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2940 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2942 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2943 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2944 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2947 Double_t sumw = w1+w2;
2951 w1 = (d2+0.5)/(d1+d2+1.);
2952 w2 = (d1+0.5)/(d1+d2+1.);
2954 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2955 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2957 // get pair of "best" hypothesys
2959 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2960 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2962 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2963 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2964 //if (track1->fFakeRatio>0) continue;
2965 RegisterClusterTracks(track1,trackID1);
2966 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2967 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2969 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2970 //if (track2->fFakeRatio>0) continue;
2972 RegisterClusterTracks(track2,trackID2);
2973 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2974 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2975 UnRegisterClusterTracks(track2,trackID2);
2977 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2978 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2979 if (nskipped>0.5) continue;
2981 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2982 if (conflict1+1<cconflict1) continue;
2983 if (conflict2+1<cconflict2) continue;
2987 for (Int_t i=0;i<6;i++){
2989 Float_t c1 =0.; // conflict coeficients
2991 if (clist1[i]&&clist2[i]){
2994 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2997 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2999 c1 = 2./TMath::Max(3.+deltan,2.);
3000 c2 = 2./TMath::Max(3.+deltan,2.);
3006 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3009 deltan = (clist1[i]->GetNz()-nz1[i]);
3011 c1 = 2./TMath::Max(3.+deltan,2.);
3018 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3021 deltan = (clist2[i]->GetNz()-nz2[i]);
3023 c2 = 2./TMath::Max(3.+deltan,2.);
3029 if (TMath::Abs(track1->GetDy(i))>0.) {
3030 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3031 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3032 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3033 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3035 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3038 if (TMath::Abs(track2->GetDy(i))>0.) {
3039 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3040 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3041 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3042 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3045 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3047 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3048 if (chi21>0) sum+=w1;
3049 if (chi22>0) sum+=w2;
3052 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3053 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3054 Double_t normchi2 = 2*conflict+sumchi2/sum;
3055 if ( normchi2 <maxchi2 ){
3058 maxconflicts = conflict;
3062 UnRegisterClusterTracks(track1,trackID1);
3065 // if (maxconflicts<4 && maxchi2<th0){
3066 if (maxchi2<th0*2.){
3067 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3068 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3069 track1->SetChi2MIP(5,maxconflicts);
3070 track1->SetChi2MIP(6,maxchi2);
3071 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3072 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3073 track1->SetChi2MIP(8,index1);
3074 fBestTrackIndex[trackID1] =index1;
3075 UpdateESDtrack(track1, AliESDtrack::kITSin);
3077 else if (track10->GetChi2MIP(0)<th1){
3078 track10->SetChi2MIP(5,maxconflicts);
3079 track10->SetChi2MIP(6,maxchi2);
3080 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3081 UpdateESDtrack(track10,AliESDtrack::kITSin);
3084 for (Int_t itrack=0;itrack<entries1;itrack++){
3085 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3086 UnRegisterClusterTracks(track,trackID1);
3089 for (Int_t itrack=0;itrack<entries2;itrack++){
3090 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3091 UnRegisterClusterTracks(track,trackID2);
3094 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3095 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3096 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3097 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3098 RegisterClusterTracks(track10,trackID1);
3100 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3101 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3102 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3103 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3104 RegisterClusterTracks(track20,trackID2);
3109 //------------------------------------------------------------------------
3110 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3111 //--------------------------------------------------------------------
3112 // This function marks clusters assigned to the track
3113 //--------------------------------------------------------------------
3114 AliTracker::UseClusters(t,from);
3116 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3117 //if (c->GetQ()>2) c->Use();
3118 if (c->GetSigmaZ2()>0.1) c->Use();
3119 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3120 //if (c->GetQ()>2) c->Use();
3121 if (c->GetSigmaZ2()>0.1) c->Use();
3124 //------------------------------------------------------------------------
3125 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3127 //------------------------------------------------------------------
3128 // add track to the list of hypothesys
3129 //------------------------------------------------------------------
3131 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3132 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3134 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3136 array = new TObjArray(10);
3137 fTrackHypothesys.AddAt(array,esdindex);
3139 array->AddLast(track);
3141 //------------------------------------------------------------------------
3142 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3144 //-------------------------------------------------------------------
3145 // compress array of track hypothesys
3146 // keep only maxsize best hypothesys
3147 //-------------------------------------------------------------------
3148 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3149 if (! (fTrackHypothesys.At(esdindex)) ) return;
3150 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3151 Int_t entries = array->GetEntriesFast();
3153 //- find preliminary besttrack as a reference
3154 Float_t minchi2=10000;
3156 AliITStrackMI * besttrack=0;
3157 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3158 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3159 if (!track) continue;
3160 Float_t chi2 = NormalizedChi2(track,0);
3162 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3163 track->SetLabel(tpcLabel);
3165 track->SetFakeRatio(1.);
3166 CookLabel(track,0.); //For comparison only
3168 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3169 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3170 if (track->GetNumberOfClusters()<maxn) continue;
3171 maxn = track->GetNumberOfClusters();
3178 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3179 delete array->RemoveAt(itrack);
3183 if (!besttrack) return;
3186 //take errors of best track as a reference
3187 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3188 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3189 for (Int_t j=0;j<6;j++) {
3190 if (besttrack->GetClIndex(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);
3198 // calculate normalized chi2
3200 Float_t * chi2 = new Float_t[entries];
3201 Int_t * index = new Int_t[entries];
3202 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3203 for (Int_t itrack=0;itrack<entries;itrack++){
3204 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3206 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3207 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3208 chi2[itrack] = track->GetChi2MIP(0);
3210 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3211 delete array->RemoveAt(itrack);
3217 TMath::Sort(entries,chi2,index,kFALSE);
3218 besttrack = (AliITStrackMI*)array->At(index[0]);
3219 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3220 for (Int_t j=0;j<6;j++){
3221 if (besttrack->GetClIndex(j)>=0){
3222 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3223 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3224 ny[j] = besttrack->GetNy(j);
3225 nz[j] = besttrack->GetNz(j);
3230 // calculate one more time with updated normalized errors
3231 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3232 for (Int_t itrack=0;itrack<entries;itrack++){
3233 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3235 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3236 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3237 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3240 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3241 delete array->RemoveAt(itrack);
3246 entries = array->GetEntriesFast();
3250 TObjArray * newarray = new TObjArray();
3251 TMath::Sort(entries,chi2,index,kFALSE);
3252 besttrack = (AliITStrackMI*)array->At(index[0]);
3255 for (Int_t j=0;j<6;j++){
3256 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3257 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3258 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3259 ny[j] = besttrack->GetNy(j);
3260 nz[j] = besttrack->GetNz(j);
3263 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3264 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3265 Float_t minn = besttrack->GetNumberOfClusters()-3;
3267 for (Int_t i=0;i<entries;i++){
3268 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3269 if (!track) continue;
3270 if (accepted>maxcut) break;
3271 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3272 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3273 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3274 delete array->RemoveAt(index[i]);
3278 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3279 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3280 if (!shortbest) accepted++;
3282 newarray->AddLast(array->RemoveAt(index[i]));
3283 for (Int_t j=0;j<6;j++){
3285 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3286 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3287 ny[j] = track->GetNy(j);
3288 nz[j] = track->GetNz(j);
3293 delete array->RemoveAt(index[i]);
3297 delete fTrackHypothesys.RemoveAt(esdindex);
3298 fTrackHypothesys.AddAt(newarray,esdindex);
3302 delete fTrackHypothesys.RemoveAt(esdindex);
3308 //------------------------------------------------------------------------
3309 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3311 //-------------------------------------------------------------
3312 // try to find best hypothesy
3313 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3314 //-------------------------------------------------------------
3315 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3316 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3317 if (!array) return 0;
3318 Int_t entries = array->GetEntriesFast();
3319 if (!entries) return 0;
3320 Float_t minchi2 = 100000;
3321 AliITStrackMI * besttrack=0;
3323 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3324 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3325 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3326 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3328 for (Int_t i=0;i<entries;i++){
3329 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3330 if (!track) continue;
3331 Float_t sigmarfi,sigmaz;
3332 GetDCASigma(track,sigmarfi,sigmaz);
3333 track->SetDnorm(0,sigmarfi);
3334 track->SetDnorm(1,sigmaz);
3336 track->SetChi2MIP(1,1000000);
3337 track->SetChi2MIP(2,1000000);
3338 track->SetChi2MIP(3,1000000);
3341 backtrack = new(backtrack) AliITStrackMI(*track);
3342 if (track->GetConstrain()) {
3343 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3344 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3345 backtrack->ResetCovariance(10.);
3347 backtrack->ResetCovariance(10.);
3349 backtrack->ResetClusters();
3351 Double_t x = original->GetX();
3352 if (!RefitAt(x,backtrack,track)) continue;
3354 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3355 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3356 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3357 track->SetChi22(GetMatchingChi2(backtrack,original));
3359 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3360 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3361 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3364 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3366 //forward track - without constraint
3367 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3368 forwardtrack->ResetClusters();
3370 RefitAt(x,forwardtrack,track);
3371 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3372 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3373 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3375 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3376 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3377 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3378 forwardtrack->SetD(0,track->GetD(0));
3379 forwardtrack->SetD(1,track->GetD(1));
3382 AliITSRecPoint* clist[6];
3383 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3384 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3387 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3388 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3389 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3390 track->SetChi2MIP(3,1000);
3393 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3395 for (Int_t ichi=0;ichi<5;ichi++){
3396 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3398 if (chi2 < minchi2){
3399 //besttrack = new AliITStrackMI(*forwardtrack);
3401 besttrack->SetLabel(track->GetLabel());
3402 besttrack->SetFakeRatio(track->GetFakeRatio());
3404 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3405 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3406 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3410 delete forwardtrack;
3412 for (Int_t i=0;i<entries;i++){
3413 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3415 if (!track) continue;
3417 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3418 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3419 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3420 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3421 delete array->RemoveAt(i);
3431 SortTrackHypothesys(esdindex,checkmax,1);
3433 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3434 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3435 besttrack = (AliITStrackMI*)array->At(0);
3436 if (!besttrack) return 0;
3437 besttrack->SetChi2MIP(8,0);
3438 fBestTrackIndex[esdindex]=0;
3439 entries = array->GetEntriesFast();
3440 AliITStrackMI *longtrack =0;
3442 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3443 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3444 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3445 if (!track->GetConstrain()) continue;
3446 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3447 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3448 if (track->GetChi2MIP(0)>4.) continue;
3449 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3452 //if (longtrack) besttrack=longtrack;
3455 AliITSRecPoint * clist[6];
3456 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3457 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3458 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3459 RegisterClusterTracks(besttrack,esdindex);
3466 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3467 if (sharedtrack>=0){
3469 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3471 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3477 if (shared>2.5) return 0;
3478 if (shared>1.0) return besttrack;
3480 // Don't sign clusters if not gold track
3482 if (!besttrack->IsGoldPrimary()) return besttrack;
3483 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3485 if (fConstraint[fPass]){
3489 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3490 for (Int_t i=0;i<6;i++){
3491 Int_t index = besttrack->GetClIndex(i);
3492 if (index<0) continue;
3493 Int_t ilayer = (index & 0xf0000000) >> 28;
3494 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3495 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3497 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3498 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3499 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3500 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3501 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3502 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3504 Bool_t cansign = kTRUE;
3505 for (Int_t itrack=0;itrack<entries; itrack++){
3506 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3507 if (!track) continue;
3508 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3509 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3515 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3516 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3517 if (!c->IsUsed()) c->Use();
3523 //------------------------------------------------------------------------
3524 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3527 // get "best" hypothesys
3530 Int_t nentries = itsTracks.GetEntriesFast();
3531 for (Int_t i=0;i<nentries;i++){
3532 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3533 if (!track) continue;
3534 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3535 if (!array) continue;
3536 if (array->GetEntriesFast()<=0) continue;
3538 AliITStrackMI* longtrack=0;
3540 Float_t maxchi2=1000;
3541 for (Int_t j=0;j<array->GetEntriesFast();j++){
3542 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3543 if (!trackHyp) continue;
3544 if (trackHyp->GetGoldV0()) {
3545 longtrack = trackHyp; //gold V0 track taken
3548 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3549 Float_t chi2 = trackHyp->GetChi2MIP(0);
3551 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3553 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3555 if (chi2 > maxchi2) continue;
3556 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3563 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3564 if (!longtrack) {longtrack = besttrack;}
3565 else besttrack= longtrack;
3569 AliITSRecPoint * clist[6];
3570 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3572 track->SetNUsed(shared);
3573 track->SetNSkipped(besttrack->GetNSkipped());
3574 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3576 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3580 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3581 //if (sharedtrack==-1) sharedtrack=0;
3582 if (sharedtrack>=0) {
3583 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3586 if (besttrack&&fAfterV0) {
3587 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3589 if (besttrack&&fConstraint[fPass])
3590 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3591 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3592 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3593 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3599 //------------------------------------------------------------------------
3600 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3601 //--------------------------------------------------------------------
3602 //This function "cooks" a track label. If label<0, this track is fake.
3603 //--------------------------------------------------------------------
3606 if ( track->GetESDtrack()) tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3608 track->SetChi2MIP(9,0);
3610 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3611 Int_t cindex = track->GetClusterIndex(i);
3612 Int_t l=(cindex & 0xf0000000) >> 28;
3613 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3615 for (Int_t ind=0;ind<3;ind++){
3617 if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3618 AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3620 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3623 Int_t nclusters = track->GetNumberOfClusters();
3624 if (nclusters > 0) //PH Some tracks don't have any cluster
3625 track->SetFakeRatio(double(nwrong)/double(nclusters));
3627 if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3629 track->SetLabel(tpcLabel);
3631 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3634 //------------------------------------------------------------------------
3635 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3638 track->SetChi2MIP(9,0);
3639 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3640 Int_t cindex = track->GetClusterIndex(i);
3641 Int_t l=(cindex & 0xf0000000) >> 28;
3642 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3643 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3645 for (Int_t ind=0;ind<3;ind++){
3646 if (cl->GetLabel(ind)==lab) isWrong=0;
3648 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3652 track->CookdEdx(low,up);
3654 //------------------------------------------------------------------------
3655 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3657 // Create some arrays
3659 if (fCoefficients) delete []fCoefficients;
3660 fCoefficients = new Float_t[ntracks*48];
3661 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3663 //------------------------------------------------------------------------
3664 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3667 // Compute predicted chi2
3669 Float_t erry,errz,covyz;
3670 Float_t theta = track->GetTgl();
3671 Float_t phi = track->GetSnp();
3672 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3673 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3674 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()));
3675 // Take into account the mis-alignment (bring track to cluster plane)
3676 Double_t xTrOrig=track->GetX();
3677 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3678 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()));
3679 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3680 // Bring the track back to detector plane in ideal geometry
3681 // [mis-alignment will be accounted for in UpdateMI()]
3682 if (!track->Propagate(xTrOrig)) return 1000.;
3684 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3685 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3687 chi2+=0.5*TMath::Min(delta/2,2.);
3688 chi2+=2.*cluster->GetDeltaProbability();
3691 track->SetNy(layer,ny);
3692 track->SetNz(layer,nz);
3693 track->SetSigmaY(layer,erry);
3694 track->SetSigmaZ(layer, errz);
3695 track->SetSigmaYZ(layer,covyz);
3696 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3697 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3701 //------------------------------------------------------------------------
3702 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3707 Int_t layer = (index & 0xf0000000) >> 28;
3708 track->SetClIndex(layer, index);
3709 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3710 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3711 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3712 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3716 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
3719 // Take into account the mis-alignment (bring track to cluster plane)
3720 Double_t xTrOrig=track->GetX();
3721 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3722 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3723 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3724 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3726 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3729 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3730 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3731 c.SetSigmaYZ(track->GetSigmaYZ(layer));
3734 Int_t updated = track->UpdateMI(&c,chi2,index);
3735 // Bring the track back to detector plane in ideal geometry
3736 if (!track->Propagate(xTrOrig)) return 0;
3738 if(!updated) AliDebug(2,"update failed");
3742 //------------------------------------------------------------------------
3743 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3746 //DCA sigmas parameterization
3747 //to be paramterized using external parameters in future
3750 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3751 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3753 //------------------------------------------------------------------------
3754 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3757 // Clusters from delta electrons?
3759 Int_t entries = clusterArray->GetEntriesFast();
3760 if (entries<4) return;
3761 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3762 Int_t layer = cluster->GetLayer();
3763 if (layer>1) return;
3765 Int_t ncandidates=0;
3766 Float_t r = (layer>0)? 7:4;
3768 for (Int_t i=0;i<entries;i++){
3769 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3770 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3771 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3772 index[ncandidates] = i; //candidate to belong to delta electron track
3774 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3775 cl0->SetDeltaProbability(1);
3781 for (Int_t i=0;i<ncandidates;i++){
3782 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3783 if (cl0->GetDeltaProbability()>0.8) continue;
3786 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3787 sumy=sumz=sumy2=sumyz=sumw=0.0;
3788 for (Int_t j=0;j<ncandidates;j++){
3790 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3792 Float_t dz = cl0->GetZ()-cl1->GetZ();
3793 Float_t dy = cl0->GetY()-cl1->GetY();
3794 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3795 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3796 y[ncl] = cl1->GetY();
3797 z[ncl] = cl1->GetZ();
3798 sumy+= y[ncl]*weight;
3799 sumz+= z[ncl]*weight;
3800 sumy2+=y[ncl]*y[ncl]*weight;
3801 sumyz+=y[ncl]*z[ncl]*weight;
3806 if (ncl<4) continue;
3807 Float_t det = sumw*sumy2 - sumy*sumy;
3809 if (TMath::Abs(det)>0.01){
3810 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3811 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3812 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3815 Float_t z0 = sumyz/sumy;
3816 delta = TMath::Abs(cl0->GetZ()-z0);
3819 cl0->SetDeltaProbability(1-20.*delta);
3823 //------------------------------------------------------------------------
3824 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3829 track->UpdateESDtrack(flags);
3830 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3831 if (oldtrack) delete oldtrack;
3832 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3833 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3834 printf("Problem\n");
3837 //------------------------------------------------------------------------
3838 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3840 // Get nearest upper layer close to the point xr.
3841 // rough approximation
3843 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3844 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3846 for (Int_t i=0;i<6;i++){
3847 if (radius<kRadiuses[i]){
3854 //------------------------------------------------------------------------
3855 void AliITStrackerMI::BuildMaterialLUT(TString material) {
3856 //--------------------------------------------------------------------
3857 // Fill a look-up table with mean material
3858 //--------------------------------------------------------------------
3862 Double_t point1[3],point2[3];
3863 Double_t phi,cosphi,sinphi,z;
3864 // 0-5 layers, 6 pipe, 7-8 shields
3865 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
3866 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
3868 Int_t ifirst=0,ilast=0;
3869 if(material.Contains("Pipe")) {
3871 } else if(material.Contains("Shields")) {
3873 } else if(material.Contains("Layers")) {
3876 Error("BuildMaterialLUT","Wrong layer name\n");
3879 for(Int_t imat=ifirst; imat<=ilast; imat++) {
3880 Double_t param[5]={0.,0.,0.,0.,0.};
3881 for (Int_t i=0; i<n; i++) {
3882 phi = 2.*TMath::Pi()*gRandom->Rndm();
3883 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
3884 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
3885 point1[0] = rmin[imat]*cosphi;
3886 point1[1] = rmin[imat]*sinphi;
3888 point2[0] = rmax[imat]*cosphi;
3889 point2[1] = rmax[imat]*sinphi;
3891 AliTracker::MeanMaterialBudget(point1,point2,mparam);
3892 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
3894 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
3896 fxOverX0Layer[imat] = param[1];
3897 fxTimesRhoLayer[imat] = param[0]*param[4];
3898 } else if(imat==6) {
3899 fxOverX0Pipe = param[1];
3900 fxTimesRhoPipe = param[0]*param[4];
3901 } else if(imat==7) {
3902 fxOverX0Shield[0] = param[1];
3903 fxTimesRhoShield[0] = param[0]*param[4];
3904 } else if(imat==8) {
3905 fxOverX0Shield[1] = param[1];
3906 fxTimesRhoShield[1] = param[0]*param[4];
3910 printf("%s\n",material.Data());
3911 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
3912 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
3913 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
3914 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
3915 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
3916 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
3917 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
3918 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
3919 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
3923 //------------------------------------------------------------------------
3924 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
3925 TString direction) {
3926 //-------------------------------------------------------------------
3927 // Propagate beyond beam pipe and correct for material
3928 // (material budget in different ways according to fUseTGeo value)
3929 // Add time if going outward (PropagateTo or PropagateToTGeo)
3930 //-------------------------------------------------------------------
3932 // Define budget mode:
3933 // 0: material from AliITSRecoParam (hard coded)
3934 // 1: material from TGeo in one step (on the fly)
3935 // 2: material from lut
3936 // 3: material from TGeo in one step (same for all hypotheses)
3949 if(fTrackingPhase.Contains("Clusters2Tracks"))
3950 { mode=3; } else { mode=1; }
3953 if(fTrackingPhase.Contains("Clusters2Tracks"))
3954 { mode=3; } else { mode=2; }
3960 if(fTrackingPhase.Contains("Default")) mode=0;
3962 Int_t index=fCurrentEsdTrack;
3964 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
3965 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
3967 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
3969 Double_t xOverX0,x0,lengthTimesMeanDensity;
3973 xOverX0 = AliITSRecoParam::GetdPipe();
3974 x0 = AliITSRecoParam::GetX0Be();
3975 lengthTimesMeanDensity = xOverX0*x0;
3976 lengthTimesMeanDensity *= dir;
3977 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
3980 if (!t->PropagateToTGeo(xToGo,1)) return 0;
3983 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
3984 xOverX0 = fxOverX0Pipe;
3985 lengthTimesMeanDensity = fxTimesRhoPipe;
3986 lengthTimesMeanDensity *= dir;
3987 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
3990 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
3991 if(fxOverX0PipeTrks[index]<0) {
3992 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
3993 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
3994 ((1.-t->GetSnp())*(1.+t->GetSnp())));
3995 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
3996 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
3999 xOverX0 = fxOverX0PipeTrks[index];
4000 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4001 lengthTimesMeanDensity *= dir;
4002 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4008 //------------------------------------------------------------------------
4009 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4011 TString direction) {
4012 //-------------------------------------------------------------------
4013 // Propagate beyond SPD or SDD shield and correct for material
4014 // (material budget in different ways according to fUseTGeo value)
4015 // Add time if going outward (PropagateTo or PropagateToTGeo)
4016 //-------------------------------------------------------------------
4018 // Define budget mode:
4019 // 0: material from AliITSRecoParam (hard coded)
4020 // 1: material from TGeo in steps of X cm (on the fly)
4021 // X = AliITSRecoParam::GetStepSizeTGeo()
4022 // 2: material from lut
4023 // 3: material from TGeo in one step (same for all hypotheses)
4036 if(fTrackingPhase.Contains("Clusters2Tracks"))
4037 { mode=3; } else { mode=1; }
4040 if(fTrackingPhase.Contains("Clusters2Tracks"))
4041 { mode=3; } else { mode=2; }
4047 if(fTrackingPhase.Contains("Default")) mode=0;
4049 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4051 Int_t shieldindex=0;
4052 if (shield.Contains("SDD")) { // SDDouter
4053 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4055 } else if (shield.Contains("SPD")) { // SPDouter
4056 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4059 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4063 // do nothing if we are already beyond the shield
4064 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4065 if(dir<0 && rTrack > rToGo) return 1; // going outward
4066 if(dir>0 && rTrack < rToGo) return 1; // going inward
4070 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4072 Int_t index=2*fCurrentEsdTrack+shieldindex;
4074 Double_t xOverX0,x0,lengthTimesMeanDensity;
4079 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4080 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4081 lengthTimesMeanDensity = xOverX0*x0;
4082 lengthTimesMeanDensity *= dir;
4083 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4086 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4087 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4090 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4091 xOverX0 = fxOverX0Shield[shieldindex];
4092 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4093 lengthTimesMeanDensity *= dir;
4094 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4097 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4098 if(fxOverX0ShieldTrks[index]<0) {
4099 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4100 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4101 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4102 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4103 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4106 xOverX0 = fxOverX0ShieldTrks[index];
4107 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4108 lengthTimesMeanDensity *= dir;
4109 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4115 //------------------------------------------------------------------------
4116 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4118 Double_t oldGlobXYZ[3],
4119 TString direction) {
4120 //-------------------------------------------------------------------
4121 // Propagate beyond layer and correct for material
4122 // (material budget in different ways according to fUseTGeo value)
4123 // Add time if going outward (PropagateTo or PropagateToTGeo)
4124 //-------------------------------------------------------------------
4126 // Define budget mode:
4127 // 0: material from AliITSRecoParam (hard coded)
4128 // 1: material from TGeo in stepsof X cm (on the fly)
4129 // X = AliITSRecoParam::GetStepSizeTGeo()
4130 // 2: material from lut
4131 // 3: material from TGeo in one step (same for all hypotheses)
4144 if(fTrackingPhase.Contains("Clusters2Tracks"))
4145 { mode=3; } else { mode=1; }
4148 if(fTrackingPhase.Contains("Clusters2Tracks"))
4149 { mode=3; } else { mode=2; }
4155 if(fTrackingPhase.Contains("Default")) mode=0;
4157 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4159 Double_t r=fgLayers[layerindex].GetR();
4160 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4162 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4164 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4166 Int_t index=6*fCurrentEsdTrack+layerindex;
4169 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4172 // back before material (no correction)
4174 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4175 if (!t->GetLocalXat(rOld,xOld)) return 0;
4176 if (!t->Propagate(xOld)) return 0;
4180 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4181 lengthTimesMeanDensity = xOverX0*x0;
4182 lengthTimesMeanDensity *= dir;
4183 // Bring the track beyond the material
4184 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4187 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4188 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4191 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4192 xOverX0 = fxOverX0Layer[layerindex];
4193 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4194 lengthTimesMeanDensity *= dir;
4195 // Bring the track beyond the material
4196 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4199 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4200 if(fxOverX0LayerTrks[index]<0) {
4201 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4202 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4203 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4204 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4205 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4206 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4209 xOverX0 = fxOverX0LayerTrks[index];
4210 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4211 lengthTimesMeanDensity *= dir;
4212 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4219 //------------------------------------------------------------------------
4220 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4221 //-----------------------------------------------------------------
4222 // Initialize LUT for storing material for each prolonged track
4223 //-----------------------------------------------------------------
4224 fxOverX0PipeTrks = new Float_t[ntracks];
4225 fxTimesRhoPipeTrks = new Float_t[ntracks];
4226 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4227 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4228 fxOverX0LayerTrks = new Float_t[ntracks*6];
4229 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4231 for(Int_t i=0; i<ntracks; i++) {
4232 fxOverX0PipeTrks[i] = -1.;
4233 fxTimesRhoPipeTrks[i] = -1.;
4235 for(Int_t j=0; j<ntracks*2; j++) {
4236 fxOverX0ShieldTrks[j] = -1.;
4237 fxTimesRhoShieldTrks[j] = -1.;
4239 for(Int_t k=0; k<ntracks*6; k++) {
4240 fxOverX0LayerTrks[k] = -1.;
4241 fxTimesRhoLayerTrks[k] = -1.;
4248 //------------------------------------------------------------------------
4249 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4250 //-----------------------------------------------------------------
4251 // Delete LUT for storing material for each prolonged track
4252 //-----------------------------------------------------------------
4253 if(fxOverX0PipeTrks) {
4254 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4256 if(fxOverX0ShieldTrks) {
4257 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4260 if(fxOverX0LayerTrks) {
4261 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4263 if(fxTimesRhoPipeTrks) {
4264 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4266 if(fxTimesRhoShieldTrks) {
4267 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4269 if(fxTimesRhoLayerTrks) {
4270 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4274 //------------------------------------------------------------------------
4275 void AliITStrackerMI::SetForceSkippingOfLayer() {
4276 //-----------------------------------------------------------------
4277 // Check if we are forced to skip layers
4278 // either we set to skip them in RecoParam
4279 // or they were off during data-taking
4280 //-----------------------------------------------------------------
4282 const AliEventInfo *eventInfo = GetEventInfo();
4284 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4285 fForceSkippingOfLayer[l] = 0;
4287 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4291 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4293 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4294 } else if(l==2 || l==3) {
4295 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4297 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4303 //------------------------------------------------------------------------
4304 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4305 Int_t ilayer,Int_t idet) const {
4306 //-----------------------------------------------------------------
4307 // This method is used to decide whether to allow a prolongation
4308 // without clusters, because we want to skip the layer.
4309 // In this case the return value is > 0:
4310 // return 1: the user requested to skip a layer
4311 // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
4312 //-----------------------------------------------------------------
4314 if (ForceSkippingOfLayer(ilayer)) return 1;
4316 if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4317 // check if track will cross SPD outer layer
4318 Double_t phiAtSPD2,zAtSPD2;
4319 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4320 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4326 //------------------------------------------------------------------------
4327 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4328 Int_t ilayer,Int_t idet,
4329 Double_t dz,Double_t dy,
4330 Bool_t noClusters) const {
4331 //-----------------------------------------------------------------
4332 // This method is used to decide whether to allow a prolongation
4333 // without clusters, because there is a dead zone in the road.
4334 // In this case the return value is > 0:
4335 // return 1: dead zone at z=0,+-7cm in SPD
4336 // return 2: all road is "bad" (dead or noisy) from the OCDB
4337 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4338 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4339 //-----------------------------------------------------------------
4341 // check dead zones at z=0,+-7cm in the SPD
4342 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4343 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4344 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4345 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4346 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4347 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4348 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4349 for (Int_t i=0; i<3; i++)
4350 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4351 AliDebug(2,Form("crack SPD %d",ilayer));
4356 // check bad zones from OCDB
4357 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4359 if (idet<0) return 0;
4361 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4364 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4365 if (ilayer==0 || ilayer==1) { // ---------- SPD
4367 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4369 detSizeFactorX *= 2.;
4370 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4373 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4374 if (detType==2) segm->SetLayer(ilayer+1);
4375 Float_t detSizeX = detSizeFactorX*segm->Dx();
4376 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4378 // check if the road overlaps with bad chips
4380 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
4381 Float_t zlocmin = zloc-dz;
4382 Float_t zlocmax = zloc+dz;
4383 Float_t xlocmin = xloc-dy;
4384 Float_t xlocmax = xloc+dy;
4385 Int_t chipsInRoad[100];
4387 // check if road goes out of detector
4388 Bool_t touchNeighbourDet=kFALSE;
4389 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
4390 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
4391 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
4392 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
4393 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));
4395 // check if this detector is bad
4397 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4398 if(!touchNeighbourDet) {
4399 return 2; // all detectors in road are bad
4401 return 3; // at least one is bad
4405 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4406 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4407 if (!nChipsInRoad) return 0;
4409 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4410 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4411 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4412 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4413 if (det.IsChipBad(chipsInRoad[iCh])) {
4421 if(!touchNeighbourDet) {
4422 AliDebug(2,"all bad in road");
4423 return 2; // all chips in road are bad
4425 return 3; // at least a bad chip in road
4430 AliDebug(2,"at least a bad in road");
4431 return 3; // at least a bad chip in road
4435 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4436 || !noClusters) return 0;
4438 // There are no clusters in road: check if there is at least
4439 // a bad SPD pixel or SDD anode or SSD strips on both sides
4441 Int_t idetInITS=idet;
4442 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4444 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4445 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4448 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4452 //------------------------------------------------------------------------
4453 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4454 const AliITStrackMI *track,
4455 Float_t &xloc,Float_t &zloc) const {
4456 //-----------------------------------------------------------------
4457 // Gives position of track in local module ref. frame
4458 //-----------------------------------------------------------------
4463 if(idet<0) return kFALSE;
4465 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4467 Int_t lad = Int_t(idet/ndet) + 1;
4469 Int_t det = idet - (lad-1)*ndet + 1;
4471 Double_t xyzGlob[3],xyzLoc[3];
4473 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4474 // take into account the misalignment: xyz at real detector plane
4475 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4477 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4479 xloc = (Float_t)xyzLoc[0];
4480 zloc = (Float_t)xyzLoc[2];
4484 //------------------------------------------------------------------------
4485 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4487 // Method to be optimized further:
4488 // Aim: decide whether a track can be used for PlaneEff evaluation
4489 // the decision is taken based on the track quality at the layer under study
4490 // no information on the clusters on this layer has to be used
4491 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4492 // the cut is done on number of sigmas from the boundaries
4494 // Input: Actual track, layer [0,5] under study
4496 // Return: kTRUE if this is a good track
4498 // it will apply a pre-selection to obtain good quality tracks.
4499 // Here also you will have the possibility to put a control on the
4500 // impact point of the track on the basic block, in order to exclude border regions
4501 // this will be done by calling a proper method of the AliITSPlaneEff class.
4503 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4504 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4506 Int_t index[AliITSgeomTGeo::kNLayers];
4508 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4510 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4511 index[k]=clusters[k];
4515 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4516 AliITSlayer &layer=fgLayers[ilayer];
4517 Double_t r=layer.GetR();
4518 AliITStrackMI tmp(*track);
4520 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4522 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
4523 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4524 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4525 if (tmp.GetClIndex(lay)>=0) ncl++;
4527 Bool_t nextout = kFALSE;
4528 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4529 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4530 Bool_t nextin = kFALSE;
4531 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4532 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4533 if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4535 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4536 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4537 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4538 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4542 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4543 Int_t idet=layer.FindDetectorIndex(phi,z);
4544 if(idet<0) { AliInfo(Form("cannot find detector"));
4547 // here check if it has good Chi Square.
4549 //propagate to the intersection with the detector plane
4550 const AliITSdetector &det=layer.GetDetector(idet);
4551 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4555 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4556 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4557 if(key>fPlaneEff->Nblock()) return kFALSE;
4558 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4559 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4561 // DEFINITION OF SEARCH ROAD FOR accepting a track
4563 //For the time being they are hard-wired, later on from AliITSRecoParam
4564 // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
4565 // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
4568 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4569 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4570 // done for RecPoints
4572 // exclude tracks at boundary between detectors
4573 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4574 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4575 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4576 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4577 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4579 if ( (locx-dx < blockXmn+boundaryWidth) ||
4580 (locx+dx > blockXmx-boundaryWidth) ||
4581 (locz-dz < blockZmn+boundaryWidth) ||
4582 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4585 //------------------------------------------------------------------------
4586 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4588 // This Method has to be optimized! For the time-being it uses the same criteria
4589 // as those used in the search of extra clusters for overlapping modules.
4591 // Method Purpose: estabilish whether a track has produced a recpoint or not
4592 // in the layer under study (For Plane efficiency)
4594 // inputs: AliITStrackMI* track (pointer to a usable track)
4596 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4597 // traversed by this very track. In details:
4598 // - if a cluster can be associated to the track then call
4599 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4600 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4603 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4604 AliITSlayer &layer=fgLayers[ilayer];
4605 Double_t r=layer.GetR();
4606 AliITStrackMI tmp(*track);
4610 if (!tmp.GetPhiZat(r,phi,z)) return;
4611 Int_t idet=layer.FindDetectorIndex(phi,z);
4613 if(idet<0) { AliInfo(Form("cannot find detector"));
4617 //propagate to the intersection with the detector plane
4618 const AliITSdetector &det=layer.GetDetector(idet);
4619 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4623 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4625 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4626 TMath::Sqrt(tmp.GetSigmaZ2() +
4627 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4628 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4629 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4630 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4631 TMath::Sqrt(tmp.GetSigmaY2() +
4632 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4633 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4634 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4636 // road in global (rphi,z) [i.e. in tracking ref. system]
4637 Double_t zmin = tmp.GetZ() - dz;
4638 Double_t zmax = tmp.GetZ() + dz;
4639 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4640 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4642 // select clusters in road
4643 layer.SelectClusters(zmin,zmax,ymin,ymax);
4645 // Define criteria for track-cluster association
4646 Double_t msz = tmp.GetSigmaZ2() +
4647 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4648 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4649 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4650 Double_t msy = tmp.GetSigmaY2() +
4651 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4652 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4653 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4654 if (tmp.GetConstrain()) {
4655 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4656 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4658 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4659 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4661 msz = 1./msz; // 1/RoadZ^2
4662 msy = 1./msy; // 1/RoadY^2
4665 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4667 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4668 //Double_t tolerance=0.2;
4669 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4670 idetc = cl->GetDetectorIndex();
4671 if(idet!=idetc) continue;
4672 //Int_t ilay = cl->GetLayer();
4674 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4675 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4677 Double_t chi2=tmp.GetPredictedChi2(cl);
4678 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4682 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4684 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4685 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4686 if(key>fPlaneEff->Nblock()) return;
4687 Bool_t found=kFALSE;
4690 while ((cl=layer.GetNextCluster(clidx))!=0) {
4691 idetc = cl->GetDetectorIndex();
4692 if(idet!=idetc) continue;
4693 // here real control to see whether the cluster can be associated to the track.
4694 // cluster not associated to track
4695 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4696 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
4697 // calculate track-clusters chi2
4698 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4699 // in particular, the error associated to the cluster
4700 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4702 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4704 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4705 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4706 // track->SetExtraModule(ilayer,idetExtra);
4708 if(!fPlaneEff->UpDatePlaneEff(found,key))
4709 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4710 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4711 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
4712 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4713 Int_t cltype[2]={-999,-999};
4717 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4718 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4721 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4722 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4723 cltype[0]=layer.GetCluster(ci)->GetNy();
4724 cltype[1]=layer.GetCluster(ci)->GetNz();
4726 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4727 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
4728 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4729 // It is computed properly by calling the method
4730 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4732 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4733 //if (tmp.PropagateTo(x,0.,0.)) {
4734 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4735 AliCluster c(*layer.GetCluster(ci));
4736 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4737 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4738 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
4739 clu[2]=TMath::Sqrt(c.GetSigmaY2());
4740 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4743 fPlaneEff->FillHistos(key,found,tr,clu,cltype);