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 "AliITSRecPointContainer.h"
54 #include "AliITSgeomTGeo.h"
55 #include "AliITSReconstructor.h"
56 #include "AliITSClusterParam.h"
57 #include "AliITSsegmentation.h"
58 #include "AliITSCalibration.h"
59 #include "AliITSPlaneEffSPD.h"
60 #include "AliITSPlaneEffSDD.h"
61 #include "AliITSPlaneEffSSD.h"
62 #include "AliITSV0Finder.h"
63 #include "AliITStrackerMI.h"
64 #include "AliMathBase.h"
66 ClassImp(AliITStrackerMI)
68 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
70 AliITStrackerMI::AliITStrackerMI():AliTracker(),
80 fLastLayerToTrackTo(0),
83 fTrackingPhase("Default"),
89 fxTimesRhoPipeTrks(0),
90 fxOverX0ShieldTrks(0),
91 fxTimesRhoShieldTrks(0),
93 fxTimesRhoLayerTrks(0),
100 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
101 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
102 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
104 //------------------------------------------------------------------------
105 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
106 fI(AliITSgeomTGeo::GetNLayers()),
115 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
118 fTrackingPhase("Default"),
124 fxTimesRhoPipeTrks(0),
125 fxOverX0ShieldTrks(0),
126 fxTimesRhoShieldTrks(0),
127 fxOverX0LayerTrks(0),
128 fxTimesRhoLayerTrks(0),
130 fITSChannelStatus(0),
133 //--------------------------------------------------------------------
134 //This is the AliITStrackerMI constructor
135 //--------------------------------------------------------------------
137 AliWarning("\"geom\" is actually a dummy argument !");
143 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
144 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
145 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
147 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
148 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
149 Double_t poff=TMath::ATan2(y,x);
151 Double_t r=TMath::Sqrt(x*x + y*y);
153 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
154 r += TMath::Sqrt(x*x + y*y);
155 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
156 r += TMath::Sqrt(x*x + y*y);
157 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
158 r += TMath::Sqrt(x*x + y*y);
161 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
163 for (Int_t j=1; j<nlad+1; j++) {
164 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
165 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
166 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
168 Double_t txyz[3]={0.};
169 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
170 m.LocalToMaster(txyz,xyz);
171 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
172 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
174 if (phi<0) phi+=TMath::TwoPi();
175 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
177 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
178 new(&det) AliITSdetector(r,phi);
179 // compute the real radius (with misalignment)
180 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
182 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
183 mmisal.LocalToMaster(txyz,xyz);
184 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
185 det.SetRmisal(rmisal);
187 } // end loop on detectors
188 } // end loop on ladders
189 fForceSkippingOfLayer[i] = 0;
190 } // end loop on layers
193 fI=AliITSgeomTGeo::GetNLayers();
196 fConstraint[0]=1; fConstraint[1]=0;
198 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
199 AliITSReconstructor::GetRecoParam()->GetYVdef(),
200 AliITSReconstructor::GetRecoParam()->GetZVdef()};
201 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
202 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
203 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
204 SetVertex(xyzVtx,ersVtx);
206 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
207 for (Int_t i=0;i<100000;i++){
208 fBestTrackIndex[i]=0;
211 // store positions of centre of SPD modules (in z)
213 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
214 fSPDdetzcentre[0] = tr[2];
215 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
216 fSPDdetzcentre[1] = tr[2];
217 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
218 fSPDdetzcentre[2] = tr[2];
219 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
220 fSPDdetzcentre[3] = tr[2];
222 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
223 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
224 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
228 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
229 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
231 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
233 // only for plane efficiency evaluation
234 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
235 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
236 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
237 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
238 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
239 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
240 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
241 else fPlaneEff = new AliITSPlaneEffSSD();
242 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
243 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
244 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
247 //------------------------------------------------------------------------
248 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
250 fBestTrack(tracker.fBestTrack),
251 fTrackToFollow(tracker.fTrackToFollow),
252 fTrackHypothesys(tracker.fTrackHypothesys),
253 fBestHypothesys(tracker.fBestHypothesys),
254 fOriginal(tracker.fOriginal),
255 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
256 fPass(tracker.fPass),
257 fAfterV0(tracker.fAfterV0),
258 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
259 fCoefficients(tracker.fCoefficients),
261 fTrackingPhase(tracker.fTrackingPhase),
262 fUseTGeo(tracker.fUseTGeo),
263 fNtracks(tracker.fNtracks),
264 fxOverX0Pipe(tracker.fxOverX0Pipe),
265 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
267 fxTimesRhoPipeTrks(0),
268 fxOverX0ShieldTrks(0),
269 fxTimesRhoShieldTrks(0),
270 fxOverX0LayerTrks(0),
271 fxTimesRhoLayerTrks(0),
272 fDebugStreamer(tracker.fDebugStreamer),
273 fITSChannelStatus(tracker.fITSChannelStatus),
274 fkDetTypeRec(tracker.fkDetTypeRec),
275 fPlaneEff(tracker.fPlaneEff) {
279 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
282 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
283 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
286 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
287 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
290 //------------------------------------------------------------------------
291 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
292 //Assignment operator
293 this->~AliITStrackerMI();
294 new(this) AliITStrackerMI(tracker);
297 //------------------------------------------------------------------------
298 AliITStrackerMI::~AliITStrackerMI()
303 if (fCoefficients) delete [] fCoefficients;
304 DeleteTrksMaterialLUT();
305 if (fDebugStreamer) {
306 //fDebugStreamer->Close();
307 delete fDebugStreamer;
309 if(fITSChannelStatus) delete fITSChannelStatus;
310 if(fPlaneEff) delete fPlaneEff;
312 //------------------------------------------------------------------------
313 void AliITStrackerMI::ReadBadFromDetTypeRec() {
314 //--------------------------------------------------------------------
315 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
317 //--------------------------------------------------------------------
319 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
321 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
323 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
326 if(fITSChannelStatus) delete fITSChannelStatus;
327 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
329 // ITS detectors and chips
330 Int_t i=0,j=0,k=0,ndet=0;
331 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
332 Int_t nBadDetsPerLayer=0;
333 ndet=AliITSgeomTGeo::GetNDetectors(i);
334 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
335 for (k=1; k<ndet+1; k++) {
336 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
337 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
338 if(det.IsBad()) {nBadDetsPerLayer++;}
339 } // end loop on detectors
340 } // end loop on ladders
341 Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
342 } // end loop on layers
346 //------------------------------------------------------------------------
347 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
348 //--------------------------------------------------------------------
349 //This function loads ITS clusters
350 //--------------------------------------------------------------------
352 TClonesArray *clusters = NULL;
353 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
354 clusters=rpcont->FetchClusters(0,cTree);
355 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
356 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
359 Int_t i=0,j=0,ndet=0;
361 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
362 ndet=fgLayers[i].GetNdetectors();
363 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
364 for (; j<jmax; j++) {
365 // if (!cTree->GetEvent(j)) continue;
366 clusters = rpcont->UncheckedGetClusters(j);
367 if(!clusters)continue;
368 Int_t ncl=clusters->GetEntriesFast();
369 SignDeltas(clusters,GetZ());
372 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
373 detector=c->GetDetectorIndex();
375 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
377 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
379 AliWarning(Form("Too many clusters on layer %d!",i));
384 // add dead zone "virtual" cluster in SPD, if there is a cluster within
385 // zwindow cm from the dead zone
386 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
387 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
388 Int_t lab[4] = {0,0,0,detector};
389 Int_t info[3] = {0,0,i};
390 Float_t q = 0.; // this identifies virtual clusters
391 Float_t hit[5] = {xdead,
393 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
394 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
396 Bool_t local = kTRUE;
397 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
398 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
399 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
400 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
401 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
402 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
403 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
404 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
405 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
406 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
407 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
408 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
409 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
410 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
411 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
412 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
413 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
414 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
415 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
417 } // "virtual" clusters in SPD
421 fgLayers[i].ResetRoad(); //road defined by the cluster density
422 fgLayers[i].SortClusters();
425 // check whether we have to skip some layers
426 SetForceSkippingOfLayer();
430 //------------------------------------------------------------------------
431 void AliITStrackerMI::UnloadClusters() {
432 //--------------------------------------------------------------------
433 //This function unloads ITS clusters
434 //--------------------------------------------------------------------
435 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
437 //------------------------------------------------------------------------
438 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
439 //--------------------------------------------------------------------
440 // Publishes all pointers to clusters known to the tracker into the
441 // passed object array.
442 // The ownership is not transfered - the caller is not expected to delete
444 //--------------------------------------------------------------------
446 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
447 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
448 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
455 //------------------------------------------------------------------------
456 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
457 //--------------------------------------------------------------------
458 // Correction for the material between the TPC and the ITS
459 //--------------------------------------------------------------------
460 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
461 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
462 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
463 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
464 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
465 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
466 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
467 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
469 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
475 //------------------------------------------------------------------------
476 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
477 //--------------------------------------------------------------------
478 // This functions reconstructs ITS tracks
479 // The clusters must be already loaded !
480 //--------------------------------------------------------------------
482 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
484 fTrackingPhase="Clusters2Tracks";
486 TObjArray itsTracks(15000);
488 fEsd = event; // store pointer to the esd
490 // temporary (for cosmics)
491 if(event->GetVertex()) {
492 TString title = event->GetVertex()->GetTitle();
493 if(title.Contains("cosmics")) {
494 Double_t xyz[3]={GetX(),GetY(),GetZ()};
495 Double_t exyz[3]={0.1,0.1,0.1};
501 {/* Read ESD tracks */
502 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
503 Int_t nentr=event->GetNumberOfTracks();
505 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
507 AliESDtrack *esd=event->GetTrack(nentr);
508 // ---- for debugging:
509 //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); }
511 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
512 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
513 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
514 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
517 t=new AliITStrackMI(*esd);
518 } catch (const Char_t *msg) {
519 //Warning("Clusters2Tracks",msg);
523 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
524 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
527 // look at the ESD mass hypothesys !
528 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
530 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
532 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
533 //track - can be V0 according to TPC
535 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
539 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
543 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
547 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
552 t->SetReconstructed(kFALSE);
553 itsTracks.AddLast(t);
554 fOriginal.AddLast(t);
556 } /* End Read ESD tracks */
560 Int_t nentr=itsTracks.GetEntriesFast();
561 fTrackHypothesys.Expand(nentr);
562 fBestHypothesys.Expand(nentr);
563 MakeCoefficients(nentr);
564 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
566 // THE TWO TRACKING PASSES
567 for (fPass=0; fPass<2; fPass++) {
568 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
569 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
570 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
571 if (t==0) continue; //this track has been already tracked
572 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
573 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
574 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
575 if (fConstraint[fPass]) {
576 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
577 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
580 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
581 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
583 ResetTrackToFollow(*t);
586 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
589 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
591 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
592 if (!besttrack) continue;
593 besttrack->SetLabel(tpcLabel);
594 // besttrack->CookdEdx();
596 besttrack->SetFakeRatio(1.);
597 CookLabel(besttrack,0.); //For comparison only
598 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
600 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
602 t->SetReconstructed(kTRUE);
604 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
606 GetBestHypothesysMIP(itsTracks);
607 } // end loop on the two tracking passes
609 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
610 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
615 Int_t entries = fTrackHypothesys.GetEntriesFast();
616 for (Int_t ientry=0; ientry<entries; ientry++) {
617 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
618 if (array) array->Delete();
619 delete fTrackHypothesys.RemoveAt(ientry);
622 fTrackHypothesys.Delete();
623 fBestHypothesys.Delete();
625 delete [] fCoefficients;
627 DeleteTrksMaterialLUT();
629 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
631 fTrackingPhase="Default";
635 //------------------------------------------------------------------------
636 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
637 //--------------------------------------------------------------------
638 // This functions propagates reconstructed ITS tracks back
639 // The clusters must be loaded !
640 //--------------------------------------------------------------------
641 fTrackingPhase="PropagateBack";
642 Int_t nentr=event->GetNumberOfTracks();
643 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
646 for (Int_t i=0; i<nentr; i++) {
647 AliESDtrack *esd=event->GetTrack(i);
649 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
650 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
654 t=new AliITStrackMI(*esd);
655 } catch (const Char_t *msg) {
656 //Warning("PropagateBack",msg);
660 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
662 ResetTrackToFollow(*t);
665 // propagate to vertex [SR, GSI 17.02.2003]
666 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
667 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
668 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
669 fTrackToFollow.StartTimeIntegral();
670 // from vertex to outside pipe
671 CorrectForPipeMaterial(&fTrackToFollow,"outward");
673 // Start time integral and add distance from current position to vertex
674 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
675 fTrackToFollow.GetXYZ(xyzTrk);
677 for (Int_t icoord=0; icoord<3; icoord++) {
678 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
681 fTrackToFollow.StartTimeIntegral();
682 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
684 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
685 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
686 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
690 fTrackToFollow.SetLabel(t->GetLabel());
691 //fTrackToFollow.CookdEdx();
692 CookLabel(&fTrackToFollow,0.); //For comparison only
693 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
694 //UseClusters(&fTrackToFollow);
700 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
702 fTrackingPhase="Default";
706 //------------------------------------------------------------------------
707 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
708 //--------------------------------------------------------------------
709 // This functions refits ITS tracks using the
710 // "inward propagated" TPC tracks
711 // The clusters must be loaded !
712 //--------------------------------------------------------------------
713 fTrackingPhase="RefitInward";
715 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
717 Int_t nentr=event->GetNumberOfTracks();
718 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
721 for (Int_t i=0; i<nentr; i++) {
722 AliESDtrack *esd=event->GetTrack(i);
724 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
725 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
726 if (esd->GetStatus()&AliESDtrack::kTPCout)
727 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
731 t=new AliITStrackMI(*esd);
732 } catch (const Char_t *msg) {
733 //Warning("RefitInward",msg);
737 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
738 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
743 ResetTrackToFollow(*t);
744 fTrackToFollow.ResetClusters();
746 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
747 fTrackToFollow.ResetCovariance(10.);
750 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
751 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
753 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
754 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
755 AliDebug(2," refit OK");
756 fTrackToFollow.SetLabel(t->GetLabel());
757 // fTrackToFollow.CookdEdx();
758 CookdEdx(&fTrackToFollow);
760 CookLabel(&fTrackToFollow,0.0); //For comparison only
763 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
764 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
765 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
766 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
767 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
768 Double_t r[3]={0.,0.,0.};
770 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
777 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
779 fTrackingPhase="Default";
783 //------------------------------------------------------------------------
784 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
785 //--------------------------------------------------------------------
786 // Return pointer to a given cluster
787 //--------------------------------------------------------------------
788 Int_t l=(index & 0xf0000000) >> 28;
789 Int_t c=(index & 0x0fffffff) >> 00;
790 return fgLayers[l].GetCluster(c);
792 //------------------------------------------------------------------------
793 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
794 //--------------------------------------------------------------------
795 // Get track space point with index i
796 //--------------------------------------------------------------------
798 Int_t l=(index & 0xf0000000) >> 28;
799 Int_t c=(index & 0x0fffffff) >> 00;
800 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
801 Int_t idet = cl->GetDetectorIndex();
805 cl->GetGlobalXYZ(xyz);
806 cl->GetGlobalCov(cov);
808 p.SetCharge(cl->GetQ());
809 p.SetDriftTime(cl->GetDriftTime());
810 p.SetChargeRatio(cl->GetChargeRatio());
811 p.SetClusterType(cl->GetClusterType());
812 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
815 iLayer = AliGeomManager::kSPD1;
818 iLayer = AliGeomManager::kSPD2;
821 iLayer = AliGeomManager::kSDD1;
824 iLayer = AliGeomManager::kSDD2;
827 iLayer = AliGeomManager::kSSD1;
830 iLayer = AliGeomManager::kSSD2;
833 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
836 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
837 p.SetVolumeID((UShort_t)volid);
840 //------------------------------------------------------------------------
841 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
842 AliTrackPoint& p, const AliESDtrack *t) {
843 //--------------------------------------------------------------------
844 // Get track space point with index i
845 // (assign error estimated during the tracking)
846 //--------------------------------------------------------------------
848 Int_t l=(index & 0xf0000000) >> 28;
849 Int_t c=(index & 0x0fffffff) >> 00;
850 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
851 Int_t idet = cl->GetDetectorIndex();
853 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
855 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
857 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
858 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
859 Double_t alpha = t->GetAlpha();
860 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
861 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
862 phi += alpha-det.GetPhi();
863 Float_t tgphi = TMath::Tan(phi);
865 Float_t tgl = t->GetTgl(); // tgl about const along track
866 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
868 Float_t errtrky,errtrkz,covyz;
869 Bool_t addMisalErr=kFALSE;
870 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
874 cl->GetGlobalXYZ(xyz);
875 // cl->GetGlobalCov(cov);
876 Float_t pos[3] = {0.,0.,0.};
877 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
878 tmpcl.GetGlobalCov(cov);
881 p.SetCharge(cl->GetQ());
882 p.SetDriftTime(cl->GetDriftTime());
883 p.SetChargeRatio(cl->GetChargeRatio());
884 p.SetClusterType(cl->GetClusterType());
886 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
889 iLayer = AliGeomManager::kSPD1;
892 iLayer = AliGeomManager::kSPD2;
895 iLayer = AliGeomManager::kSDD1;
898 iLayer = AliGeomManager::kSDD2;
901 iLayer = AliGeomManager::kSSD1;
904 iLayer = AliGeomManager::kSSD2;
907 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
910 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
912 p.SetVolumeID((UShort_t)volid);
915 //------------------------------------------------------------------------
916 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
918 //--------------------------------------------------------------------
919 // Follow prolongation tree
920 //--------------------------------------------------------------------
922 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
923 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
926 AliESDtrack * esd = otrack->GetESDtrack();
927 if (esd->GetV0Index(0)>0) {
928 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
929 // mapping of ESD track is different as ITS track in Containers
930 // Need something more stable
931 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
932 for (Int_t i=0;i<3;i++){
933 Int_t index = esd->GetV0Index(i);
935 AliESDv0 * vertex = fEsd->GetV0(index);
936 if (vertex->GetStatus()<0) continue; // rejected V0
938 if (esd->GetSign()>0) {
939 vertex->SetIndex(0,esdindex);
941 vertex->SetIndex(1,esdindex);
945 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
947 bestarray = new TObjArray(5);
948 fBestHypothesys.AddAt(bestarray,esdindex);
952 //setup tree of the prolongations
954 static AliITStrackMI tracks[7][100];
955 AliITStrackMI *currenttrack;
956 static AliITStrackMI currenttrack1;
957 static AliITStrackMI currenttrack2;
958 static AliITStrackMI backuptrack;
960 Int_t nindexes[7][100];
961 Float_t normalizedchi2[100];
962 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
963 otrack->SetNSkipped(0);
964 new (&(tracks[6][0])) AliITStrackMI(*otrack);
966 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
967 Int_t modstatus = 1; // found
971 // follow prolongations
972 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
973 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
976 AliITSlayer &layer=fgLayers[ilayer];
977 Double_t r = layer.GetR();
983 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
985 if (ntracks[ilayer]>=100) break;
986 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
987 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
988 if (ntracks[ilayer]>15+ilayer){
989 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
990 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
993 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
995 // material between SSD and SDD, SDD and SPD
997 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
999 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1003 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1004 Int_t idet=layer.FindDetectorIndex(phi,z);
1006 Double_t trackGlobXYZ1[3];
1007 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1009 // Get the budget to the primary vertex for the current track being prolonged
1010 Double_t budgetToPrimVertex = GetEffectiveThickness();
1012 // check if we allow a prolongation without point
1013 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1015 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1016 // propagate to the layer radius
1017 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1018 if(!vtrack->Propagate(xToGo)) continue;
1019 // apply correction for material of the current layer
1020 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1021 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1022 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1023 vtrack->SetClIndex(ilayer,-1);
1024 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1025 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1026 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1028 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1033 // track outside layer acceptance in z
1034 if (idet<0) continue;
1036 //propagate to the intersection with the detector plane
1037 const AliITSdetector &det=layer.GetDetector(idet);
1038 new(¤ttrack2) AliITStrackMI(currenttrack1);
1039 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1040 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1041 currenttrack1.SetDetectorIndex(idet);
1042 currenttrack2.SetDetectorIndex(idet);
1043 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1046 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1048 // road in global (rphi,z) [i.e. in tracking ref. system]
1049 Double_t zmin,zmax,ymin,ymax;
1050 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1052 // select clusters in road
1053 layer.SelectClusters(zmin,zmax,ymin,ymax);
1054 //********************
1056 // Define criteria for track-cluster association
1057 Double_t msz = currenttrack1.GetSigmaZ2() +
1058 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1059 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1060 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1061 Double_t msy = currenttrack1.GetSigmaY2() +
1062 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1063 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1064 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1066 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1067 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1069 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1070 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1072 msz = 1./msz; // 1/RoadZ^2
1073 msy = 1./msy; // 1/RoadY^2
1077 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1079 const AliITSRecPoint *cl=0;
1081 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1082 Bool_t deadzoneSPD=kFALSE;
1083 currenttrack = ¤ttrack1;
1085 // check if the road contains a dead zone
1086 Bool_t noClusters = kFALSE;
1087 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1088 if (noClusters) AliDebug(2,"no clusters in road");
1089 Double_t dz=0.5*(zmax-zmin);
1090 Double_t dy=0.5*(ymax-ymin);
1091 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1092 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1093 // create a prolongation without clusters (check also if there are no clusters in the road)
1096 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1097 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1098 updatetrack->SetClIndex(ilayer,-1);
1100 modstatus = 5; // no cls in road
1101 } else if (dead==1) {
1102 modstatus = 7; // holes in z in SPD
1103 } else if (dead==2 || dead==3 || dead==4) {
1104 modstatus = 2; // dead from OCDB
1106 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1107 // apply correction for material of the current layer
1108 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1109 if (constrain) { // apply vertex constrain
1110 updatetrack->SetConstrain(constrain);
1111 Bool_t isPrim = kTRUE;
1112 if (ilayer<4) { // check that it's close to the vertex
1113 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1114 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1115 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1116 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1117 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1119 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1121 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1123 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1124 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1126 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1127 updatetrack->SetDeadZoneProbability(ilayer,1.);
1128 } else if (dead==4) { // at least a single dead channel from OCDB
1129 updatetrack->SetDeadZoneProbability(ilayer,0.);
1136 // loop over clusters in the road
1137 while ((cl=layer.GetNextCluster(clidx))!=0) {
1138 if (ntracks[ilayer]>95) break; //space for skipped clusters
1139 Bool_t changedet =kFALSE;
1140 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1141 Int_t idetc=cl->GetDetectorIndex();
1143 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1144 // take into account misalignment (bring track to real detector plane)
1145 Double_t xTrOrig = currenttrack->GetX();
1146 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1147 // a first cut on track-cluster distance
1148 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1149 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1150 { // cluster not associated to track
1151 AliDebug(2,"not associated");
1154 // bring track back to ideal detector plane
1155 if (!currenttrack->Propagate(xTrOrig)) continue;
1156 } else { // have to move track to cluster's detector
1157 const AliITSdetector &detc=layer.GetDetector(idetc);
1158 // a first cut on track-cluster distance
1160 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1161 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1162 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1163 continue; // cluster not associated to track
1165 new (&backuptrack) AliITStrackMI(currenttrack2);
1167 currenttrack =¤ttrack2;
1168 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1169 new (currenttrack) AliITStrackMI(backuptrack);
1173 currenttrack->SetDetectorIndex(idetc);
1174 // Get again the budget to the primary vertex
1175 // for the current track being prolonged, if had to change detector
1176 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1179 // calculate track-clusters chi2
1180 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1182 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1183 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1184 if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1185 if (ntracks[ilayer]>=100) continue;
1186 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1187 updatetrack->SetClIndex(ilayer,-1);
1188 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1190 if (cl->GetQ()!=0) { // real cluster
1191 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1192 AliDebug(2,"update failed");
1195 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1196 modstatus = 1; // found
1197 } else { // virtual cluster in dead zone
1198 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1199 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1200 modstatus = 7; // holes in z in SPD
1204 Float_t xlocnewdet,zlocnewdet;
1205 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1206 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1209 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1211 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1213 // apply correction for material of the current layer
1214 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1216 if (constrain) { // apply vertex constrain
1217 updatetrack->SetConstrain(constrain);
1218 Bool_t isPrim = kTRUE;
1219 if (ilayer<4) { // check that it's close to the vertex
1220 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1221 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1222 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1223 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1224 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1226 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1227 } //apply vertex constrain
1229 } // create new hypothesis
1231 AliDebug(2,"chi2 too large");
1234 } // loop over possible prolongations
1236 // allow one prolongation without clusters
1237 if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1238 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1239 // apply correction for material of the current layer
1240 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1241 vtrack->SetClIndex(ilayer,-1);
1242 modstatus = 3; // skipped
1243 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1244 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1245 vtrack->IncrementNSkipped();
1250 } // loop over tracks in layer ilayer+1
1252 //loop over track candidates for the current layer
1258 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1259 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1260 if (normalizedchi2[itrack] <
1261 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1265 if (constrain) { // constrain
1266 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1268 } else { // !constrain
1269 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1274 // sort tracks by increasing normalized chi2
1275 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1276 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1277 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1278 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1279 } // end loop over layers
1283 // Now select tracks to be kept
1285 Int_t max = constrain ? 20 : 5;
1287 // tracks that reach layer 0 (SPD inner)
1288 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1289 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1290 if (track.GetNumberOfClusters()<2) continue;
1291 if (!constrain && track.GetNormChi2(0) >
1292 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1295 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1298 // tracks that reach layer 1 (SPD outer)
1299 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1300 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1301 if (track.GetNumberOfClusters()<4) continue;
1302 if (!constrain && track.GetNormChi2(1) >
1303 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1304 if (constrain) track.IncrementNSkipped();
1306 track.SetD(0,track.GetD(GetX(),GetY()));
1307 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1308 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1309 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1312 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1315 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1317 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1318 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1319 if (track.GetNumberOfClusters()<3) continue;
1320 if (!constrain && track.GetNormChi2(2) >
1321 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1322 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1324 track.SetD(0,track.GetD(GetX(),GetY()));
1325 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1326 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1327 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1330 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1336 // register best track of each layer - important for V0 finder
1338 for (Int_t ilayer=0;ilayer<5;ilayer++){
1339 if (ntracks[ilayer]==0) continue;
1340 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1341 if (track.GetNumberOfClusters()<1) continue;
1342 CookLabel(&track,0);
1343 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1347 // update TPC V0 information
1349 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1350 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1351 for (Int_t i=0;i<3;i++){
1352 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1353 if (index==0) break;
1354 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1355 if (vertex->GetStatus()<0) continue; // rejected V0
1357 if (otrack->GetSign()>0) {
1358 vertex->SetIndex(0,esdindex);
1361 vertex->SetIndex(1,esdindex);
1363 //find nearest layer with track info
1364 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1365 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1366 Int_t nearest = nearestold;
1367 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1368 if (ntracks[nearest]==0){
1373 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1374 if (nearestold<5&&nearest<5){
1375 Bool_t accept = track.GetNormChi2(nearest)<10;
1377 if (track.GetSign()>0) {
1378 vertex->SetParamP(track);
1379 vertex->Update(fprimvertex);
1380 //vertex->SetIndex(0,track.fESDtrack->GetID());
1381 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1383 vertex->SetParamN(track);
1384 vertex->Update(fprimvertex);
1385 //vertex->SetIndex(1,track.fESDtrack->GetID());
1386 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1388 vertex->SetStatus(vertex->GetStatus()+1);
1390 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1397 //------------------------------------------------------------------------
1398 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1400 //--------------------------------------------------------------------
1403 return fgLayers[layer];
1405 //------------------------------------------------------------------------
1406 AliITStrackerMI::AliITSlayer::AliITSlayer():
1436 //--------------------------------------------------------------------
1437 //default AliITSlayer constructor
1438 //--------------------------------------------------------------------
1439 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1440 fClusterWeight[i]=0;
1441 fClusterTracks[0][i]=-1;
1442 fClusterTracks[1][i]=-1;
1443 fClusterTracks[2][i]=-1;
1444 fClusterTracks[3][i]=-1;
1447 //------------------------------------------------------------------------
1448 AliITStrackerMI::AliITSlayer::
1449 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1478 //--------------------------------------------------------------------
1479 //main AliITSlayer constructor
1480 //--------------------------------------------------------------------
1481 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1482 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1484 //------------------------------------------------------------------------
1485 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1487 fPhiOffset(layer.fPhiOffset),
1488 fNladders(layer.fNladders),
1489 fZOffset(layer.fZOffset),
1490 fNdetectors(layer.fNdetectors),
1491 fDetectors(layer.fDetectors),
1496 fClustersCs(layer.fClustersCs),
1497 fClusterIndexCs(layer.fClusterIndexCs),
1501 fCurrentSlice(layer.fCurrentSlice),
1509 fAccepted(layer.fAccepted),
1511 fMaxSigmaClY(layer.fMaxSigmaClY),
1512 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1513 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1517 //------------------------------------------------------------------------
1518 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1519 //--------------------------------------------------------------------
1520 // AliITSlayer destructor
1521 //--------------------------------------------------------------------
1522 delete [] fDetectors;
1523 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1524 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1525 fClusterWeight[i]=0;
1526 fClusterTracks[0][i]=-1;
1527 fClusterTracks[1][i]=-1;
1528 fClusterTracks[2][i]=-1;
1529 fClusterTracks[3][i]=-1;
1532 //------------------------------------------------------------------------
1533 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1534 //--------------------------------------------------------------------
1535 // This function removes loaded clusters
1536 //--------------------------------------------------------------------
1537 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1538 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1539 fClusterWeight[i]=0;
1540 fClusterTracks[0][i]=-1;
1541 fClusterTracks[1][i]=-1;
1542 fClusterTracks[2][i]=-1;
1543 fClusterTracks[3][i]=-1;
1549 //------------------------------------------------------------------------
1550 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1551 //--------------------------------------------------------------------
1552 // This function reset weights of the clusters
1553 //--------------------------------------------------------------------
1554 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1555 fClusterWeight[i]=0;
1556 fClusterTracks[0][i]=-1;
1557 fClusterTracks[1][i]=-1;
1558 fClusterTracks[2][i]=-1;
1559 fClusterTracks[3][i]=-1;
1561 for (Int_t i=0; i<fN;i++) {
1562 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1563 if (cl&&cl->IsUsed()) cl->Use();
1567 //------------------------------------------------------------------------
1568 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1569 //--------------------------------------------------------------------
1570 // This function calculates the road defined by the cluster density
1571 //--------------------------------------------------------------------
1573 for (Int_t i=0; i<fN; i++) {
1574 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1576 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1578 //------------------------------------------------------------------------
1579 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1580 //--------------------------------------------------------------------
1581 //This function adds a cluster to this layer
1582 //--------------------------------------------------------------------
1583 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1589 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1591 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1592 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1593 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1594 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1595 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1596 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1599 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1600 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1601 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1602 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1606 //------------------------------------------------------------------------
1607 void AliITStrackerMI::AliITSlayer::SortClusters()
1612 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1613 Float_t *z = new Float_t[fN];
1614 Int_t * index = new Int_t[fN];
1616 fMaxSigmaClY=0.; //AD
1617 fMaxSigmaClZ=0.; //AD
1619 for (Int_t i=0;i<fN;i++){
1620 z[i] = fClusters[i]->GetZ();
1621 // save largest errors in y and z for this layer
1622 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1623 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1625 TMath::Sort(fN,z,index,kFALSE);
1626 for (Int_t i=0;i<fN;i++){
1627 clusters[i] = fClusters[index[i]];
1630 for (Int_t i=0;i<fN;i++){
1631 fClusters[i] = clusters[i];
1632 fZ[i] = fClusters[i]->GetZ();
1633 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1634 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1635 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1645 for (Int_t i=0;i<fN;i++){
1646 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1647 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1648 fClusterIndex[i] = i;
1652 fDy5 = (fYB[1]-fYB[0])/5.;
1653 fDy10 = (fYB[1]-fYB[0])/10.;
1654 fDy20 = (fYB[1]-fYB[0])/20.;
1655 for (Int_t i=0;i<6;i++) fN5[i] =0;
1656 for (Int_t i=0;i<11;i++) fN10[i]=0;
1657 for (Int_t i=0;i<21;i++) fN20[i]=0;
1659 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;}
1660 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;}
1661 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;}
1664 for (Int_t i=0;i<fN;i++)
1665 for (Int_t irot=-1;irot<=1;irot++){
1666 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1668 for (Int_t slice=0; slice<6;slice++){
1669 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1670 fClusters5[slice][fN5[slice]] = fClusters[i];
1671 fY5[slice][fN5[slice]] = curY;
1672 fZ5[slice][fN5[slice]] = fZ[i];
1673 fClusterIndex5[slice][fN5[slice]]=i;
1678 for (Int_t slice=0; slice<11;slice++){
1679 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1680 fClusters10[slice][fN10[slice]] = fClusters[i];
1681 fY10[slice][fN10[slice]] = curY;
1682 fZ10[slice][fN10[slice]] = fZ[i];
1683 fClusterIndex10[slice][fN10[slice]]=i;
1688 for (Int_t slice=0; slice<21;slice++){
1689 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1690 fClusters20[slice][fN20[slice]] = fClusters[i];
1691 fY20[slice][fN20[slice]] = curY;
1692 fZ20[slice][fN20[slice]] = fZ[i];
1693 fClusterIndex20[slice][fN20[slice]]=i;
1700 // consistency check
1702 for (Int_t i=0;i<fN-1;i++){
1708 for (Int_t slice=0;slice<21;slice++)
1709 for (Int_t i=0;i<fN20[slice]-1;i++){
1710 if (fZ20[slice][i]>fZ20[slice][i+1]){
1717 //------------------------------------------------------------------------
1718 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1719 //--------------------------------------------------------------------
1720 // This function returns the index of the nearest cluster
1721 //--------------------------------------------------------------------
1724 if (fCurrentSlice<0) {
1733 if (ncl==0) return 0;
1734 Int_t b=0, e=ncl-1, m=(b+e)/2;
1735 for (; b<e; m=(b+e)/2) {
1736 // if (z > fClusters[m]->GetZ()) b=m+1;
1737 if (z > zcl[m]) b=m+1;
1742 //------------------------------------------------------------------------
1743 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 {
1744 //--------------------------------------------------------------------
1745 // This function computes the rectangular road for this track
1746 //--------------------------------------------------------------------
1749 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1750 // take into account the misalignment: propagate track to misaligned detector plane
1751 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1753 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1754 TMath::Sqrt(track->GetSigmaZ2() +
1755 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1756 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1757 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1758 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1759 TMath::Sqrt(track->GetSigmaY2() +
1760 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1761 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1762 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1764 // track at boundary between detectors, enlarge road
1765 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1766 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1767 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1768 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1769 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1770 Float_t tgl = TMath::Abs(track->GetTgl());
1771 if (tgl > 1.) tgl=1.;
1772 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1773 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1774 Float_t snp = TMath::Abs(track->GetSnp());
1775 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1776 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1779 // add to the road a term (up to 2-3 mm) to deal with misalignments
1780 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1781 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1783 Double_t r = fgLayers[ilayer].GetR();
1784 zmin = track->GetZ() - dz;
1785 zmax = track->GetZ() + dz;
1786 ymin = track->GetY() + r*det.GetPhi() - dy;
1787 ymax = track->GetY() + r*det.GetPhi() + dy;
1789 // bring track back to idead detector plane
1790 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1794 //------------------------------------------------------------------------
1795 void AliITStrackerMI::AliITSlayer::
1796 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1797 //--------------------------------------------------------------------
1798 // This function sets the "window"
1799 //--------------------------------------------------------------------
1801 Double_t circle=2*TMath::Pi()*fR;
1807 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1808 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1809 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1810 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1811 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1813 Float_t ymiddle = (fYmin+fYmax)*0.5;
1814 if (ymiddle<fYB[0]) {
1815 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1816 } else if (ymiddle>fYB[1]) {
1817 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1823 fClustersCs = fClusters;
1824 fClusterIndexCs = fClusterIndex;
1830 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1831 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1832 if (slice<0) slice=0;
1833 if (slice>20) slice=20;
1834 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1836 fCurrentSlice=slice;
1837 fClustersCs = fClusters20[fCurrentSlice];
1838 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1839 fYcs = fY20[fCurrentSlice];
1840 fZcs = fZ20[fCurrentSlice];
1841 fNcs = fN20[fCurrentSlice];
1846 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1847 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1848 if (slice<0) slice=0;
1849 if (slice>10) slice=10;
1850 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1852 fCurrentSlice=slice;
1853 fClustersCs = fClusters10[fCurrentSlice];
1854 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1855 fYcs = fY10[fCurrentSlice];
1856 fZcs = fZ10[fCurrentSlice];
1857 fNcs = fN10[fCurrentSlice];
1862 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1863 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1864 if (slice<0) slice=0;
1865 if (slice>5) slice=5;
1866 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1868 fCurrentSlice=slice;
1869 fClustersCs = fClusters5[fCurrentSlice];
1870 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1871 fYcs = fY5[fCurrentSlice];
1872 fZcs = fZ5[fCurrentSlice];
1873 fNcs = fN5[fCurrentSlice];
1877 fI = FindClusterIndex(fZmin);
1878 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
1884 //------------------------------------------------------------------------
1885 Int_t AliITStrackerMI::AliITSlayer::
1886 FindDetectorIndex(Double_t phi, Double_t z) const {
1887 //--------------------------------------------------------------------
1888 //This function finds the detector crossed by the track
1889 //--------------------------------------------------------------------
1891 if (fZOffset<0) // old geometry
1892 dphi = -(phi-fPhiOffset);
1893 else // new geometry
1894 dphi = phi-fPhiOffset;
1897 if (dphi < 0) dphi += 2*TMath::Pi();
1898 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1899 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1900 if (np>=fNladders) np-=fNladders;
1901 if (np<0) np+=fNladders;
1904 Double_t dz=fZOffset-z;
1905 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1906 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1907 if (nz>=fNdetectors) return -1;
1908 if (nz<0) return -1;
1910 // ad hoc correction for 3rd ladder of SDD inner layer,
1911 // which is reversed (rotated by pi around local y)
1912 // this correction is OK only from AliITSv11Hybrid onwards
1913 if (GetR()>12. && GetR()<20.) { // SDD inner
1914 if(np==2) { // 3rd ladder
1915 Double_t posMod252[3];
1916 AliITSgeomTGeo::GetTranslation(252,posMod252);
1917 // check the Z coordinate of Mod 252: if negative
1918 // (old SDD geometry in AliITSv11Hybrid)
1919 // the swap of numeration whould be applied
1920 if(posMod252[2]<0.){
1921 nz = (fNdetectors-1) - nz;
1925 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1928 return np*fNdetectors + nz;
1930 //------------------------------------------------------------------------
1931 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1933 //--------------------------------------------------------------------
1934 // This function returns clusters within the "window"
1935 //--------------------------------------------------------------------
1937 if (fCurrentSlice<0) {
1938 Double_t rpi2 = 2.*fR*TMath::Pi();
1939 for (Int_t i=fI; i<fImax; i++) {
1942 if (fYmax<y) y -= rpi2;
1943 if (fYmin>y) y += rpi2;
1944 if (y<fYmin) continue;
1945 if (y>fYmax) continue;
1947 // skip clusters that are in "extended" road but they
1948 // 3sigma error does not touch the original road
1949 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
1950 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
1952 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1955 return fClusters[i];
1958 for (Int_t i=fI; i<fImax; i++) {
1959 if (fYcs[i]<fYmin) continue;
1960 if (fYcs[i]>fYmax) continue;
1961 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1962 ci=fClusterIndexCs[i];
1964 return fClustersCs[i];
1969 //------------------------------------------------------------------------
1970 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1972 //--------------------------------------------------------------------
1973 // This function returns the layer thickness at this point (units X0)
1974 //--------------------------------------------------------------------
1976 x0=AliITSRecoParam::GetX0Air();
1977 if (43<fR&&fR<45) { //SSD2
1980 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1981 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1982 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1983 for (Int_t i=0; i<12; i++) {
1984 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1985 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1989 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1990 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1994 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1995 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1998 if (37<fR&&fR<41) { //SSD1
2001 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2002 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2003 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2004 for (Int_t i=0; i<11; i++) {
2005 if (TMath::Abs(z-3.9*i)<0.15) {
2006 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2010 if (TMath::Abs(z+3.9*i)<0.15) {
2011 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2015 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2016 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2019 if (13<fR&&fR<26) { //SDD
2022 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2024 if (TMath::Abs(y-1.80)<0.55) {
2026 for (Int_t j=0; j<20; j++) {
2027 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2028 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2031 if (TMath::Abs(y+1.80)<0.55) {
2033 for (Int_t j=0; j<20; j++) {
2034 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2035 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2039 for (Int_t i=0; i<4; i++) {
2040 if (TMath::Abs(z-7.3*i)<0.60) {
2042 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2045 if (TMath::Abs(z+7.3*i)<0.60) {
2047 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2052 if (6<fR&&fR<8) { //SPD2
2053 Double_t dd=0.0063; x0=21.5;
2055 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2056 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2058 if (3<fR&&fR<5) { //SPD1
2059 Double_t dd=0.0063; x0=21.5;
2061 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2062 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2067 //------------------------------------------------------------------------
2068 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2070 fRmisal(det.fRmisal),
2072 fSinPhi(det.fSinPhi),
2073 fCosPhi(det.fCosPhi),
2079 fNChips(det.fNChips),
2080 fChipIsBad(det.fChipIsBad)
2084 //------------------------------------------------------------------------
2085 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2086 const AliITSDetTypeRec *detTypeRec)
2088 //--------------------------------------------------------------------
2089 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2090 //--------------------------------------------------------------------
2092 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2093 // while in the tracker they start from 0 for each layer
2094 for(Int_t il=0; il<ilayer; il++)
2095 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2098 if (ilayer==0 || ilayer==1) { // ---------- SPD
2100 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2102 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2105 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2109 // Get calibration from AliITSDetTypeRec
2110 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2111 calib->SetModuleIndex(idet);
2112 AliITSCalibration *calibSPDdead = 0;
2113 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2114 if (calib->IsBad() ||
2115 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2118 // printf("lay %d bad %d\n",ilayer,idet);
2121 // Get segmentation from AliITSDetTypeRec
2122 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2124 // Read info about bad chips
2125 fNChips = segm->GetMaximumChipIndex()+1;
2126 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2127 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2128 fChipIsBad = new Bool_t[fNChips];
2129 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2130 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2131 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2132 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2137 //------------------------------------------------------------------------
2138 Double_t AliITStrackerMI::GetEffectiveThickness()
2140 //--------------------------------------------------------------------
2141 // Returns the thickness between the current layer and the vertex (units X0)
2142 //--------------------------------------------------------------------
2145 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2146 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2147 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2151 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2152 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2156 Double_t xn=fgLayers[fI].GetR();
2157 for (Int_t i=0; i<fI; i++) {
2158 Double_t xi=fgLayers[i].GetR();
2159 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2165 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2166 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2169 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2170 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2174 //------------------------------------------------------------------------
2175 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2176 //-------------------------------------------------------------------
2177 // This function returns number of clusters within the "window"
2178 //--------------------------------------------------------------------
2180 for (Int_t i=fI; i<fN; i++) {
2181 const AliITSRecPoint *c=fClusters[i];
2182 if (c->GetZ() > fZmax) break;
2183 if (c->IsUsed()) continue;
2184 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2185 Double_t y=fR*det.GetPhi() + c->GetY();
2187 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2188 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2190 if (y<fYmin) continue;
2191 if (y>fYmax) continue;
2196 //------------------------------------------------------------------------
2197 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2198 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2200 //--------------------------------------------------------------------
2201 // This function refits the track "track" at the position "x" using
2202 // the clusters from "clusters"
2203 // If "extra"==kTRUE,
2204 // the clusters from overlapped modules get attached to "track"
2205 // If "planeff"==kTRUE,
2206 // special approach for plane efficiency evaluation is applyed
2207 //--------------------------------------------------------------------
2209 Int_t index[AliITSgeomTGeo::kNLayers];
2211 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2212 Int_t nc=clusters->GetNumberOfClusters();
2213 for (k=0; k<nc; k++) {
2214 Int_t idx=clusters->GetClusterIndex(k);
2215 Int_t ilayer=(idx&0xf0000000)>>28;
2219 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2221 //------------------------------------------------------------------------
2222 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2223 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2225 //--------------------------------------------------------------------
2226 // This function refits the track "track" at the position "x" using
2227 // the clusters from array
2228 // If "extra"==kTRUE,
2229 // the clusters from overlapped modules get attached to "track"
2230 // If "planeff"==kTRUE,
2231 // special approach for plane efficiency evaluation is applyed
2232 //--------------------------------------------------------------------
2233 Int_t index[AliITSgeomTGeo::kNLayers];
2235 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2237 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2238 index[k]=clusters[k];
2241 // special for cosmics: check which the innermost layer crossed
2243 Int_t innermostlayer=5;
2244 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2245 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2246 if(drphi < fgLayers[innermostlayer].GetR()) break;
2248 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2250 Int_t modstatus=1; // found
2252 Int_t from, to, step;
2253 if (xx > track->GetX()) {
2254 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2257 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2260 TString dir = (step>0 ? "outward" : "inward");
2262 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2263 AliITSlayer &layer=fgLayers[ilayer];
2264 Double_t r=layer.GetR();
2266 if (step<0 && xx>r) break;
2268 // material between SSD and SDD, SDD and SPD
2269 Double_t hI=ilayer-0.5*step;
2270 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2271 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2272 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2273 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2276 Double_t oldGlobXYZ[3];
2277 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2279 // continue if we are already beyond this layer
2280 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2281 if(step>0 && oldGlobR > r) continue; // going outward
2282 if(step<0 && oldGlobR < r) continue; // going inward
2285 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2287 Int_t idet=layer.FindDetectorIndex(phi,z);
2289 // check if we allow a prolongation without point for large-eta tracks
2290 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2292 modstatus = 4; // out in z
2293 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2294 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2297 // apply correction for material of the current layer
2298 // add time if going outward
2299 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2303 if (idet<0) return kFALSE;
2306 const AliITSdetector &det=layer.GetDetector(idet);
2307 // only for ITS-SA tracks refit
2308 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2310 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2312 track->SetDetectorIndex(idet);
2313 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2315 Double_t dz,zmin,zmax,dy,ymin,ymax;
2317 const AliITSRecPoint *clAcc=0;
2318 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2320 Int_t idx=index[ilayer];
2321 if (idx>=0) { // cluster in this layer
2322 modstatus = 6; // no refit
2323 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2325 if (idet != cl->GetDetectorIndex()) {
2326 idet=cl->GetDetectorIndex();
2327 const AliITSdetector &detc=layer.GetDetector(idet);
2328 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2329 track->SetDetectorIndex(idet);
2330 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2332 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2333 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2337 modstatus = 1; // found
2342 } else { // no cluster in this layer
2344 modstatus = 3; // skipped
2345 // Plane Eff determination:
2346 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2347 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2348 UseTrackForPlaneEff(track,ilayer);
2351 modstatus = 5; // no cls in road
2353 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2354 dz = 0.5*(zmax-zmin);
2355 dy = 0.5*(ymax-ymin);
2356 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2357 if (dead==1) modstatus = 7; // holes in z in SPD
2358 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2363 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2364 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2366 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2369 if (extra) { // search for extra clusters in overlapped modules
2370 AliITStrackV2 tmp(*track);
2371 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2372 layer.SelectClusters(zmin,zmax,ymin,ymax);
2374 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2376 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2377 Double_t tolerance=0.1;
2378 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2379 // only clusters in another module! (overlaps)
2380 idetExtra = clExtra->GetDetectorIndex();
2381 if (idet == idetExtra) continue;
2383 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2385 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2386 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2387 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2388 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2390 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2391 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2394 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2395 track->SetExtraModule(ilayer,idetExtra);
2397 } // end search for extra clusters in overlapped modules
2399 // Correct for material of the current layer
2401 // add time if going outward
2402 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2403 track->SetCheckInvariant(kTRUE);
2404 } // end loop on layers
2406 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2410 //------------------------------------------------------------------------
2411 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2414 // calculate normalized chi2
2415 // return NormalizedChi2(track,0);
2418 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2419 // track->fdEdxMismatch=0;
2420 Float_t dedxmismatch =0;
2421 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2423 for (Int_t i = 0;i<6;i++){
2424 if (track->GetClIndex(i)>=0){
2425 Float_t cerry, cerrz;
2426 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2428 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2431 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2432 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2433 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2435 cchi2+=(0.5-ratio)*10.;
2436 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2437 dedxmismatch+=(0.5-ratio)*10.;
2441 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2442 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2443 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2444 if (i<2) chi2+=2*cl->GetDeltaProbability();
2450 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2451 track->SetdEdxMismatch(dedxmismatch);
2455 for (Int_t i = 0;i<4;i++){
2456 if (track->GetClIndex(i)>=0){
2457 Float_t cerry, cerrz;
2458 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2459 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2462 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2463 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2467 for (Int_t i = 4;i<6;i++){
2468 if (track->GetClIndex(i)>=0){
2469 Float_t cerry, cerrz;
2470 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2471 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2474 Float_t cerryb, cerrzb;
2475 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2476 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2479 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2480 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2485 if (track->GetESDtrack()->GetTPCsignal()>85){
2486 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2488 chi2+=(0.5-ratio)*5.;
2491 chi2+=(ratio-2.0)*3;
2495 Double_t match = TMath::Sqrt(track->GetChi22());
2496 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2497 if (!track->GetConstrain()) {
2498 if (track->GetNumberOfClusters()>2) {
2499 match/=track->GetNumberOfClusters()-2.;
2504 if (match<0) match=0;
2506 // penalty factor for missing points (NDeadZone>0), but no penalty
2507 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2508 // or there is a dead from OCDB)
2509 Float_t deadzonefactor = 0.;
2510 if (track->GetNDeadZone()>0.) {
2511 Float_t sumDeadZoneProbability=0;
2512 for(Int_t ilay=0;ilay<6;ilay++) sumDeadZoneProbability+=track->GetDeadZoneProbability(ilay);
2513 Float_t nDeadZoneWithProbNot1=(Float_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2514 if(nDeadZoneWithProbNot1>0.) {
2515 Float_t deadZoneProbability = sumDeadZoneProbability - (Float_t)((Int_t)sumDeadZoneProbability);
2516 deadZoneProbability /= nDeadZoneWithProbNot1;
2517 deadzonefactor = 3.*(1.1-deadZoneProbability);
2521 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2522 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2523 1./(1.+track->GetNSkipped()));
2527 //------------------------------------------------------------------------
2528 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2531 // return matching chi2 between two tracks
2532 Double_t largeChi2=1000.;
2534 AliITStrackMI track3(*track2);
2535 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2537 vec(0,0)=track1->GetY() - track3.GetY();
2538 vec(1,0)=track1->GetZ() - track3.GetZ();
2539 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2540 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2541 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2544 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2545 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2546 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2547 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2548 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2550 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2551 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2552 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2553 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2555 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2556 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2557 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2559 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2560 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2562 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2565 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2566 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2569 //------------------------------------------------------------------------
2570 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2573 // return probability that given point (characterized by z position and error)
2574 // is in SPD dead zone
2576 Double_t probability = 0.;
2577 Double_t absz = TMath::Abs(zpos);
2578 Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2579 0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2580 if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2581 Double_t zmin, zmax;
2582 if (zpos<-6.) { // dead zone at z = -7
2583 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2584 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2585 } else if (zpos>6.) { // dead zone at z = +7
2586 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2587 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2588 } else if (absz<2.) { // dead zone at z = 0
2589 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2590 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2595 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2597 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2598 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2601 //------------------------------------------------------------------------
2602 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2605 // calculate normalized chi2
2607 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2609 for (Int_t i = 0;i<6;i++){
2610 if (TMath::Abs(track->GetDy(i))>0){
2611 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2612 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2615 else{chi2[i]=10000;}
2618 TMath::Sort(6,chi2,index,kFALSE);
2619 Float_t max = float(ncl)*fac-1.;
2620 Float_t sumchi=0, sumweight=0;
2621 for (Int_t i=0;i<max+1;i++){
2622 Float_t weight = (i<max)?1.:(max+1.-i);
2623 sumchi+=weight*chi2[index[i]];
2626 Double_t normchi2 = sumchi/sumweight;
2629 //------------------------------------------------------------------------
2630 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2633 // calculate normalized chi2
2634 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2637 for (Int_t i=0;i<6;i++){
2638 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2639 Double_t sy1 = forwardtrack->GetSigmaY(i);
2640 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2641 Double_t sy2 = backtrack->GetSigmaY(i);
2642 Double_t sz2 = backtrack->GetSigmaZ(i);
2643 if (i<2){ sy2=1000.;sz2=1000;}
2645 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2646 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2648 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2649 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2651 res+= nz0*nz0+ny0*ny0;
2654 if (npoints>1) return
2655 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2656 //2*forwardtrack->fNUsed+
2657 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2658 1./(1.+forwardtrack->GetNSkipped()));
2661 //------------------------------------------------------------------------
2662 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2663 //--------------------------------------------------------------------
2664 // Return pointer to a given cluster
2665 //--------------------------------------------------------------------
2666 Int_t l=(index & 0xf0000000) >> 28;
2667 Int_t c=(index & 0x0fffffff) >> 00;
2668 return fgLayers[l].GetWeight(c);
2670 //------------------------------------------------------------------------
2671 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2673 //---------------------------------------------
2674 // register track to the list
2676 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2679 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2680 Int_t index = track->GetClusterIndex(icluster);
2681 Int_t l=(index & 0xf0000000) >> 28;
2682 Int_t c=(index & 0x0fffffff) >> 00;
2683 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2684 for (Int_t itrack=0;itrack<4;itrack++){
2685 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2686 fgLayers[l].SetClusterTracks(itrack,c,id);
2692 //------------------------------------------------------------------------
2693 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2695 //---------------------------------------------
2696 // unregister track from the list
2697 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2698 Int_t index = track->GetClusterIndex(icluster);
2699 Int_t l=(index & 0xf0000000) >> 28;
2700 Int_t c=(index & 0x0fffffff) >> 00;
2701 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2702 for (Int_t itrack=0;itrack<4;itrack++){
2703 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2704 fgLayers[l].SetClusterTracks(itrack,c,-1);
2709 //------------------------------------------------------------------------
2710 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2712 //-------------------------------------------------------------
2713 //get number of shared clusters
2714 //-------------------------------------------------------------
2716 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2717 // mean number of clusters
2718 Float_t *ny = GetNy(id), *nz = GetNz(id);
2721 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2722 Int_t index = track->GetClusterIndex(icluster);
2723 Int_t l=(index & 0xf0000000) >> 28;
2724 Int_t c=(index & 0x0fffffff) >> 00;
2725 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2727 printf("problem\n");
2729 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2733 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2734 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2735 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2737 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2740 deltan = (cl->GetNz()-nz[l]);
2742 if (deltan>2.0) continue; // extended - highly probable shared cluster
2743 weight = 2./TMath::Max(3.+deltan,2.);
2745 for (Int_t itrack=0;itrack<4;itrack++){
2746 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2748 clist[l] = (AliITSRecPoint*)GetCluster(index);
2754 track->SetNUsed(shared);
2757 //------------------------------------------------------------------------
2758 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2761 // find first shared track
2763 // mean number of clusters
2764 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2766 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2767 Int_t sharedtrack=100000;
2768 Int_t tracks[24],trackindex=0;
2769 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2771 for (Int_t icluster=0;icluster<6;icluster++){
2772 if (clusterlist[icluster]<0) continue;
2773 Int_t index = clusterlist[icluster];
2774 Int_t l=(index & 0xf0000000) >> 28;
2775 Int_t c=(index & 0x0fffffff) >> 00;
2777 printf("problem\n");
2779 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2780 //if (l>3) continue;
2781 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2784 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2785 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2786 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2788 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2791 deltan = (cl->GetNz()-nz[l]);
2793 if (deltan>2.0) continue; // extended - highly probable shared cluster
2795 for (Int_t itrack=3;itrack>=0;itrack--){
2796 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2797 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2798 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2803 if (trackindex==0) return -1;
2805 sharedtrack = tracks[0];
2807 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2810 Int_t tracks2[24], cluster[24];
2811 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2814 for (Int_t i=0;i<trackindex;i++){
2815 if (tracks[i]<0) continue;
2816 tracks2[index] = tracks[i];
2818 for (Int_t j=i+1;j<trackindex;j++){
2819 if (tracks[j]<0) continue;
2820 if (tracks[j]==tracks[i]){
2828 for (Int_t i=0;i<index;i++){
2829 if (cluster[index]>max) {
2830 sharedtrack=tracks2[index];
2837 if (sharedtrack>=100000) return -1;
2839 // make list of overlaps
2841 for (Int_t icluster=0;icluster<6;icluster++){
2842 if (clusterlist[icluster]<0) continue;
2843 Int_t index = clusterlist[icluster];
2844 Int_t l=(index & 0xf0000000) >> 28;
2845 Int_t c=(index & 0x0fffffff) >> 00;
2846 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2847 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2849 if (cl->GetNy()>2) continue;
2850 if (cl->GetNz()>2) continue;
2853 if (cl->GetNy()>3) continue;
2854 if (cl->GetNz()>3) continue;
2857 for (Int_t itrack=3;itrack>=0;itrack--){
2858 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2859 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2867 //------------------------------------------------------------------------
2868 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2870 // try to find track hypothesys without conflicts
2871 // with minimal chi2;
2872 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2873 Int_t entries1 = arr1->GetEntriesFast();
2874 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2875 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2876 Int_t entries2 = arr2->GetEntriesFast();
2877 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2879 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2880 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2881 if (track10->Pt()>0.5+track20->Pt()) return track10;
2883 for (Int_t itrack=0;itrack<entries1;itrack++){
2884 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2885 UnRegisterClusterTracks(track,trackID1);
2888 for (Int_t itrack=0;itrack<entries2;itrack++){
2889 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2890 UnRegisterClusterTracks(track,trackID2);
2894 Float_t maxconflicts=6;
2895 Double_t maxchi2 =1000.;
2897 // get weight of hypothesys - using best hypothesy
2900 Int_t list1[6],list2[6];
2901 AliITSRecPoint *clist1[6], *clist2[6] ;
2902 RegisterClusterTracks(track10,trackID1);
2903 RegisterClusterTracks(track20,trackID2);
2904 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2905 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2906 UnRegisterClusterTracks(track10,trackID1);
2907 UnRegisterClusterTracks(track20,trackID2);
2910 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2911 Float_t nerry[6],nerrz[6];
2912 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2913 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2914 for (Int_t i=0;i<6;i++){
2915 if ( (erry1[i]>0) && (erry2[i]>0)) {
2916 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2917 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2919 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2920 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2922 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2923 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2924 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2927 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2928 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2929 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2937 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2938 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2939 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2940 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2942 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2943 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2944 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2946 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2947 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2948 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2951 Double_t sumw = w1+w2;
2955 w1 = (d2+0.5)/(d1+d2+1.);
2956 w2 = (d1+0.5)/(d1+d2+1.);
2958 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2959 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2961 // get pair of "best" hypothesys
2963 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2964 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2966 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2967 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2968 //if (track1->fFakeRatio>0) continue;
2969 RegisterClusterTracks(track1,trackID1);
2970 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2971 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2973 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2974 //if (track2->fFakeRatio>0) continue;
2976 RegisterClusterTracks(track2,trackID2);
2977 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2978 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2979 UnRegisterClusterTracks(track2,trackID2);
2981 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2982 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2983 if (nskipped>0.5) continue;
2985 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2986 if (conflict1+1<cconflict1) continue;
2987 if (conflict2+1<cconflict2) continue;
2991 for (Int_t i=0;i<6;i++){
2993 Float_t c1 =0.; // conflict coeficients
2995 if (clist1[i]&&clist2[i]){
2998 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3001 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3003 c1 = 2./TMath::Max(3.+deltan,2.);
3004 c2 = 2./TMath::Max(3.+deltan,2.);
3010 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3013 deltan = (clist1[i]->GetNz()-nz1[i]);
3015 c1 = 2./TMath::Max(3.+deltan,2.);
3022 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3025 deltan = (clist2[i]->GetNz()-nz2[i]);
3027 c2 = 2./TMath::Max(3.+deltan,2.);
3033 if (TMath::Abs(track1->GetDy(i))>0.) {
3034 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3035 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3036 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3037 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3039 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3042 if (TMath::Abs(track2->GetDy(i))>0.) {
3043 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3044 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3045 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3046 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3049 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3051 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3052 if (chi21>0) sum+=w1;
3053 if (chi22>0) sum+=w2;
3056 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3057 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3058 Double_t normchi2 = 2*conflict+sumchi2/sum;
3059 if ( normchi2 <maxchi2 ){
3062 maxconflicts = conflict;
3066 UnRegisterClusterTracks(track1,trackID1);
3069 // if (maxconflicts<4 && maxchi2<th0){
3070 if (maxchi2<th0*2.){
3071 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3072 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3073 track1->SetChi2MIP(5,maxconflicts);
3074 track1->SetChi2MIP(6,maxchi2);
3075 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3076 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3077 track1->SetChi2MIP(8,index1);
3078 fBestTrackIndex[trackID1] =index1;
3079 UpdateESDtrack(track1, AliESDtrack::kITSin);
3081 else if (track10->GetChi2MIP(0)<th1){
3082 track10->SetChi2MIP(5,maxconflicts);
3083 track10->SetChi2MIP(6,maxchi2);
3084 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3085 UpdateESDtrack(track10,AliESDtrack::kITSin);
3088 for (Int_t itrack=0;itrack<entries1;itrack++){
3089 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3090 UnRegisterClusterTracks(track,trackID1);
3093 for (Int_t itrack=0;itrack<entries2;itrack++){
3094 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3095 UnRegisterClusterTracks(track,trackID2);
3098 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3099 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3100 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3101 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3102 RegisterClusterTracks(track10,trackID1);
3104 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3105 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3106 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3107 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3108 RegisterClusterTracks(track20,trackID2);
3113 //------------------------------------------------------------------------
3114 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3115 //--------------------------------------------------------------------
3116 // This function marks clusters assigned to the track
3117 //--------------------------------------------------------------------
3118 AliTracker::UseClusters(t,from);
3120 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3121 //if (c->GetQ()>2) c->Use();
3122 if (c->GetSigmaZ2()>0.1) c->Use();
3123 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3124 //if (c->GetQ()>2) c->Use();
3125 if (c->GetSigmaZ2()>0.1) c->Use();
3128 //------------------------------------------------------------------------
3129 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3131 //------------------------------------------------------------------
3132 // add track to the list of hypothesys
3133 //------------------------------------------------------------------
3135 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3136 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3138 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3140 array = new TObjArray(10);
3141 fTrackHypothesys.AddAt(array,esdindex);
3143 array->AddLast(track);
3145 //------------------------------------------------------------------------
3146 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3148 //-------------------------------------------------------------------
3149 // compress array of track hypothesys
3150 // keep only maxsize best hypothesys
3151 //-------------------------------------------------------------------
3152 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3153 if (! (fTrackHypothesys.At(esdindex)) ) return;
3154 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3155 Int_t entries = array->GetEntriesFast();
3157 //- find preliminary besttrack as a reference
3158 Float_t minchi2=10000;
3160 AliITStrackMI * besttrack=0;
3161 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3162 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3163 if (!track) continue;
3164 Float_t chi2 = NormalizedChi2(track,0);
3166 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3167 track->SetLabel(tpcLabel);
3169 track->SetFakeRatio(1.);
3170 CookLabel(track,0.); //For comparison only
3172 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3173 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3174 if (track->GetNumberOfClusters()<maxn) continue;
3175 maxn = track->GetNumberOfClusters();
3182 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3183 delete array->RemoveAt(itrack);
3187 if (!besttrack) return;
3190 //take errors of best track as a reference
3191 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3192 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3193 for (Int_t j=0;j<6;j++) {
3194 if (besttrack->GetClIndex(j)>=0){
3195 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3196 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3197 ny[j] = besttrack->GetNy(j);
3198 nz[j] = besttrack->GetNz(j);
3202 // calculate normalized chi2
3204 Float_t * chi2 = new Float_t[entries];
3205 Int_t * index = new Int_t[entries];
3206 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3207 for (Int_t itrack=0;itrack<entries;itrack++){
3208 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3210 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3211 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3212 chi2[itrack] = track->GetChi2MIP(0);
3214 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3215 delete array->RemoveAt(itrack);
3221 TMath::Sort(entries,chi2,index,kFALSE);
3222 besttrack = (AliITStrackMI*)array->At(index[0]);
3223 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3224 for (Int_t j=0;j<6;j++){
3225 if (besttrack->GetClIndex(j)>=0){
3226 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3227 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3228 ny[j] = besttrack->GetNy(j);
3229 nz[j] = besttrack->GetNz(j);
3234 // calculate one more time with updated normalized errors
3235 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3236 for (Int_t itrack=0;itrack<entries;itrack++){
3237 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3239 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3240 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3241 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3244 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3245 delete array->RemoveAt(itrack);
3250 entries = array->GetEntriesFast();
3254 TObjArray * newarray = new TObjArray();
3255 TMath::Sort(entries,chi2,index,kFALSE);
3256 besttrack = (AliITStrackMI*)array->At(index[0]);
3259 for (Int_t j=0;j<6;j++){
3260 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3261 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3262 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3263 ny[j] = besttrack->GetNy(j);
3264 nz[j] = besttrack->GetNz(j);
3267 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3268 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3269 Float_t minn = besttrack->GetNumberOfClusters()-3;
3271 for (Int_t i=0;i<entries;i++){
3272 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3273 if (!track) continue;
3274 if (accepted>maxcut) break;
3275 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3276 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3277 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3278 delete array->RemoveAt(index[i]);
3282 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3283 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3284 if (!shortbest) accepted++;
3286 newarray->AddLast(array->RemoveAt(index[i]));
3287 for (Int_t j=0;j<6;j++){
3289 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3290 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3291 ny[j] = track->GetNy(j);
3292 nz[j] = track->GetNz(j);
3297 delete array->RemoveAt(index[i]);
3301 delete fTrackHypothesys.RemoveAt(esdindex);
3302 fTrackHypothesys.AddAt(newarray,esdindex);
3306 delete fTrackHypothesys.RemoveAt(esdindex);
3312 //------------------------------------------------------------------------
3313 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3315 //-------------------------------------------------------------
3316 // try to find best hypothesy
3317 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3318 //-------------------------------------------------------------
3319 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3320 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3321 if (!array) return 0;
3322 Int_t entries = array->GetEntriesFast();
3323 if (!entries) return 0;
3324 Float_t minchi2 = 100000;
3325 AliITStrackMI * besttrack=0;
3327 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3328 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3329 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3330 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3332 for (Int_t i=0;i<entries;i++){
3333 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3334 if (!track) continue;
3335 Float_t sigmarfi,sigmaz;
3336 GetDCASigma(track,sigmarfi,sigmaz);
3337 track->SetDnorm(0,sigmarfi);
3338 track->SetDnorm(1,sigmaz);
3340 track->SetChi2MIP(1,1000000);
3341 track->SetChi2MIP(2,1000000);
3342 track->SetChi2MIP(3,1000000);
3345 backtrack = new(backtrack) AliITStrackMI(*track);
3346 if (track->GetConstrain()) {
3347 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3348 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3349 backtrack->ResetCovariance(10.);
3351 backtrack->ResetCovariance(10.);
3353 backtrack->ResetClusters();
3355 Double_t x = original->GetX();
3356 if (!RefitAt(x,backtrack,track)) continue;
3358 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3359 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3360 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3361 track->SetChi22(GetMatchingChi2(backtrack,original));
3363 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3364 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3365 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3368 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3370 //forward track - without constraint
3371 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3372 forwardtrack->ResetClusters();
3374 RefitAt(x,forwardtrack,track);
3375 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3376 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3377 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3379 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3380 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3381 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3382 forwardtrack->SetD(0,track->GetD(0));
3383 forwardtrack->SetD(1,track->GetD(1));
3386 AliITSRecPoint* clist[6];
3387 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3388 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3391 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3392 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3393 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3394 track->SetChi2MIP(3,1000);
3397 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3399 for (Int_t ichi=0;ichi<5;ichi++){
3400 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3402 if (chi2 < minchi2){
3403 //besttrack = new AliITStrackMI(*forwardtrack);
3405 besttrack->SetLabel(track->GetLabel());
3406 besttrack->SetFakeRatio(track->GetFakeRatio());
3408 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3409 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3410 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3414 delete forwardtrack;
3416 for (Int_t i=0;i<entries;i++){
3417 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3419 if (!track) continue;
3421 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3422 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3423 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3424 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3425 delete array->RemoveAt(i);
3435 SortTrackHypothesys(esdindex,checkmax,1);
3437 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3438 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3439 besttrack = (AliITStrackMI*)array->At(0);
3440 if (!besttrack) return 0;
3441 besttrack->SetChi2MIP(8,0);
3442 fBestTrackIndex[esdindex]=0;
3443 entries = array->GetEntriesFast();
3444 AliITStrackMI *longtrack =0;
3446 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3447 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3448 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3449 if (!track->GetConstrain()) continue;
3450 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3451 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3452 if (track->GetChi2MIP(0)>4.) continue;
3453 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3456 //if (longtrack) besttrack=longtrack;
3459 AliITSRecPoint * clist[6];
3460 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3461 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3462 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3463 RegisterClusterTracks(besttrack,esdindex);
3470 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3471 if (sharedtrack>=0){
3473 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3475 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3481 if (shared>2.5) return 0;
3482 if (shared>1.0) return besttrack;
3484 // Don't sign clusters if not gold track
3486 if (!besttrack->IsGoldPrimary()) return besttrack;
3487 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3489 if (fConstraint[fPass]){
3493 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3494 for (Int_t i=0;i<6;i++){
3495 Int_t index = besttrack->GetClIndex(i);
3496 if (index<0) continue;
3497 Int_t ilayer = (index & 0xf0000000) >> 28;
3498 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3499 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3501 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3502 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3503 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3504 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3505 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3506 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3508 Bool_t cansign = kTRUE;
3509 for (Int_t itrack=0;itrack<entries; itrack++){
3510 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3511 if (!track) continue;
3512 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3513 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3519 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3520 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3521 if (!c->IsUsed()) c->Use();
3527 //------------------------------------------------------------------------
3528 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3531 // get "best" hypothesys
3534 Int_t nentries = itsTracks.GetEntriesFast();
3535 for (Int_t i=0;i<nentries;i++){
3536 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3537 if (!track) continue;
3538 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3539 if (!array) continue;
3540 if (array->GetEntriesFast()<=0) continue;
3542 AliITStrackMI* longtrack=0;
3544 Float_t maxchi2=1000;
3545 for (Int_t j=0;j<array->GetEntriesFast();j++){
3546 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3547 if (!trackHyp) continue;
3548 if (trackHyp->GetGoldV0()) {
3549 longtrack = trackHyp; //gold V0 track taken
3552 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3553 Float_t chi2 = trackHyp->GetChi2MIP(0);
3555 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3557 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3559 if (chi2 > maxchi2) continue;
3560 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3567 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3568 if (!longtrack) {longtrack = besttrack;}
3569 else besttrack= longtrack;
3573 AliITSRecPoint * clist[6];
3574 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3576 track->SetNUsed(shared);
3577 track->SetNSkipped(besttrack->GetNSkipped());
3578 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3580 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3584 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3585 //if (sharedtrack==-1) sharedtrack=0;
3586 if (sharedtrack>=0) {
3587 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3590 if (besttrack&&fAfterV0) {
3591 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3593 if (besttrack&&fConstraint[fPass])
3594 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3595 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3596 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3597 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3603 //------------------------------------------------------------------------
3604 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3605 //--------------------------------------------------------------------
3606 //This function "cooks" a track label. If label<0, this track is fake.
3607 //--------------------------------------------------------------------
3610 if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
3612 track->SetChi2MIP(9,0);
3614 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3615 Int_t cindex = track->GetClusterIndex(i);
3616 Int_t l=(cindex & 0xf0000000) >> 28;
3617 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3619 for (Int_t ind=0;ind<3;ind++){
3620 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3621 AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3623 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3626 Int_t nclusters = track->GetNumberOfClusters();
3627 if (nclusters > 0) //PH Some tracks don't have any cluster
3628 track->SetFakeRatio(double(nwrong)/double(nclusters));
3629 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3630 track->SetLabel(-tpcLabel);
3632 track->SetLabel(tpcLabel);
3634 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3637 //------------------------------------------------------------------------
3638 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3641 track->SetChi2MIP(9,0);
3642 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3643 Int_t cindex = track->GetClusterIndex(i);
3644 Int_t l=(cindex & 0xf0000000) >> 28;
3645 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3646 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3648 for (Int_t ind=0;ind<3;ind++){
3649 if (cl->GetLabel(ind)==lab) isWrong=0;
3651 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3655 track->CookdEdx(low,up);
3657 //------------------------------------------------------------------------
3658 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3660 // Create some arrays
3662 if (fCoefficients) delete []fCoefficients;
3663 fCoefficients = new Float_t[ntracks*48];
3664 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3666 //------------------------------------------------------------------------
3667 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3670 // Compute predicted chi2
3672 Float_t erry,errz,covyz;
3673 Float_t theta = track->GetTgl();
3674 Float_t phi = track->GetSnp();
3675 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3676 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3677 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()));
3678 // Take into account the mis-alignment (bring track to cluster plane)
3679 Double_t xTrOrig=track->GetX();
3680 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3681 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()));
3682 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3683 // Bring the track back to detector plane in ideal geometry
3684 // [mis-alignment will be accounted for in UpdateMI()]
3685 if (!track->Propagate(xTrOrig)) return 1000.;
3687 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3688 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3690 chi2+=0.5*TMath::Min(delta/2,2.);
3691 chi2+=2.*cluster->GetDeltaProbability();
3694 track->SetNy(layer,ny);
3695 track->SetNz(layer,nz);
3696 track->SetSigmaY(layer,erry);
3697 track->SetSigmaZ(layer, errz);
3698 track->SetSigmaYZ(layer,covyz);
3699 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3700 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3704 //------------------------------------------------------------------------
3705 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3710 Int_t layer = (index & 0xf0000000) >> 28;
3711 track->SetClIndex(layer, index);
3712 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3713 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3714 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3715 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3719 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
3722 // Take into account the mis-alignment (bring track to cluster plane)
3723 Double_t xTrOrig=track->GetX();
3724 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3725 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3726 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3727 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3729 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3732 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3733 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3734 c.SetSigmaYZ(track->GetSigmaYZ(layer));
3737 Int_t updated = track->UpdateMI(&c,chi2,index);
3738 // Bring the track back to detector plane in ideal geometry
3739 if (!track->Propagate(xTrOrig)) return 0;
3741 if(!updated) AliDebug(2,"update failed");
3745 //------------------------------------------------------------------------
3746 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3749 //DCA sigmas parameterization
3750 //to be paramterized using external parameters in future
3753 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3754 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3756 //------------------------------------------------------------------------
3757 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3760 // Clusters from delta electrons?
3762 Int_t entries = clusterArray->GetEntriesFast();
3763 if (entries<4) return;
3764 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3765 Int_t layer = cluster->GetLayer();
3766 if (layer>1) return;
3768 Int_t ncandidates=0;
3769 Float_t r = (layer>0)? 7:4;
3771 for (Int_t i=0;i<entries;i++){
3772 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3773 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3774 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3775 index[ncandidates] = i; //candidate to belong to delta electron track
3777 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3778 cl0->SetDeltaProbability(1);
3784 for (Int_t i=0;i<ncandidates;i++){
3785 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3786 if (cl0->GetDeltaProbability()>0.8) continue;
3789 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3790 sumy=sumz=sumy2=sumyz=sumw=0.0;
3791 for (Int_t j=0;j<ncandidates;j++){
3793 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3795 Float_t dz = cl0->GetZ()-cl1->GetZ();
3796 Float_t dy = cl0->GetY()-cl1->GetY();
3797 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3798 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3799 y[ncl] = cl1->GetY();
3800 z[ncl] = cl1->GetZ();
3801 sumy+= y[ncl]*weight;
3802 sumz+= z[ncl]*weight;
3803 sumy2+=y[ncl]*y[ncl]*weight;
3804 sumyz+=y[ncl]*z[ncl]*weight;
3809 if (ncl<4) continue;
3810 Float_t det = sumw*sumy2 - sumy*sumy;
3812 if (TMath::Abs(det)>0.01){
3813 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3814 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3815 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3818 Float_t z0 = sumyz/sumy;
3819 delta = TMath::Abs(cl0->GetZ()-z0);
3822 cl0->SetDeltaProbability(1-20.*delta);
3826 //------------------------------------------------------------------------
3827 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3832 track->UpdateESDtrack(flags);
3833 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3834 if (oldtrack) delete oldtrack;
3835 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3836 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3837 printf("Problem\n");
3840 //------------------------------------------------------------------------
3841 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3843 // Get nearest upper layer close to the point xr.
3844 // rough approximation
3846 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3847 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3849 for (Int_t i=0;i<6;i++){
3850 if (radius<kRadiuses[i]){
3857 //------------------------------------------------------------------------
3858 void AliITStrackerMI::BuildMaterialLUT(TString material) {
3859 //--------------------------------------------------------------------
3860 // Fill a look-up table with mean material
3861 //--------------------------------------------------------------------
3865 Double_t point1[3],point2[3];
3866 Double_t phi,cosphi,sinphi,z;
3867 // 0-5 layers, 6 pipe, 7-8 shields
3868 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
3869 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
3871 Int_t ifirst=0,ilast=0;
3872 if(material.Contains("Pipe")) {
3874 } else if(material.Contains("Shields")) {
3876 } else if(material.Contains("Layers")) {
3879 Error("BuildMaterialLUT","Wrong layer name\n");
3882 for(Int_t imat=ifirst; imat<=ilast; imat++) {
3883 Double_t param[5]={0.,0.,0.,0.,0.};
3884 for (Int_t i=0; i<n; i++) {
3885 phi = 2.*TMath::Pi()*gRandom->Rndm();
3886 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
3887 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
3888 point1[0] = rmin[imat]*cosphi;
3889 point1[1] = rmin[imat]*sinphi;
3891 point2[0] = rmax[imat]*cosphi;
3892 point2[1] = rmax[imat]*sinphi;
3894 AliTracker::MeanMaterialBudget(point1,point2,mparam);
3895 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
3897 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
3899 fxOverX0Layer[imat] = param[1];
3900 fxTimesRhoLayer[imat] = param[0]*param[4];
3901 } else if(imat==6) {
3902 fxOverX0Pipe = param[1];
3903 fxTimesRhoPipe = param[0]*param[4];
3904 } else if(imat==7) {
3905 fxOverX0Shield[0] = param[1];
3906 fxTimesRhoShield[0] = param[0]*param[4];
3907 } else if(imat==8) {
3908 fxOverX0Shield[1] = param[1];
3909 fxTimesRhoShield[1] = param[0]*param[4];
3913 printf("%s\n",material.Data());
3914 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
3915 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
3916 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
3917 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
3918 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
3919 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
3920 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
3921 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
3922 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
3926 //------------------------------------------------------------------------
3927 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
3928 TString direction) {
3929 //-------------------------------------------------------------------
3930 // Propagate beyond beam pipe and correct for material
3931 // (material budget in different ways according to fUseTGeo value)
3932 // Add time if going outward (PropagateTo or PropagateToTGeo)
3933 //-------------------------------------------------------------------
3935 // Define budget mode:
3936 // 0: material from AliITSRecoParam (hard coded)
3937 // 1: material from TGeo in one step (on the fly)
3938 // 2: material from lut
3939 // 3: material from TGeo in one step (same for all hypotheses)
3952 if(fTrackingPhase.Contains("Clusters2Tracks"))
3953 { mode=3; } else { mode=1; }
3956 if(fTrackingPhase.Contains("Clusters2Tracks"))
3957 { mode=3; } else { mode=2; }
3963 if(fTrackingPhase.Contains("Default")) mode=0;
3965 Int_t index=fCurrentEsdTrack;
3967 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
3968 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
3970 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
3972 Double_t xOverX0,x0,lengthTimesMeanDensity;
3976 xOverX0 = AliITSRecoParam::GetdPipe();
3977 x0 = AliITSRecoParam::GetX0Be();
3978 lengthTimesMeanDensity = xOverX0*x0;
3979 lengthTimesMeanDensity *= dir;
3980 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
3983 if (!t->PropagateToTGeo(xToGo,1)) return 0;
3986 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
3987 xOverX0 = fxOverX0Pipe;
3988 lengthTimesMeanDensity = fxTimesRhoPipe;
3989 lengthTimesMeanDensity *= dir;
3990 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
3993 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
3994 if(fxOverX0PipeTrks[index]<0) {
3995 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
3996 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
3997 ((1.-t->GetSnp())*(1.+t->GetSnp())));
3998 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
3999 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4002 xOverX0 = fxOverX0PipeTrks[index];
4003 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4004 lengthTimesMeanDensity *= dir;
4005 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4011 //------------------------------------------------------------------------
4012 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4014 TString direction) {
4015 //-------------------------------------------------------------------
4016 // Propagate beyond SPD or SDD shield and correct for material
4017 // (material budget in different ways according to fUseTGeo value)
4018 // Add time if going outward (PropagateTo or PropagateToTGeo)
4019 //-------------------------------------------------------------------
4021 // Define budget mode:
4022 // 0: material from AliITSRecoParam (hard coded)
4023 // 1: material from TGeo in steps of X cm (on the fly)
4024 // X = AliITSRecoParam::GetStepSizeTGeo()
4025 // 2: material from lut
4026 // 3: material from TGeo in one step (same for all hypotheses)
4039 if(fTrackingPhase.Contains("Clusters2Tracks"))
4040 { mode=3; } else { mode=1; }
4043 if(fTrackingPhase.Contains("Clusters2Tracks"))
4044 { mode=3; } else { mode=2; }
4050 if(fTrackingPhase.Contains("Default")) mode=0;
4052 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4054 Int_t shieldindex=0;
4055 if (shield.Contains("SDD")) { // SDDouter
4056 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4058 } else if (shield.Contains("SPD")) { // SPDouter
4059 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4062 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4066 // do nothing if we are already beyond the shield
4067 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4068 if(dir<0 && rTrack > rToGo) return 1; // going outward
4069 if(dir>0 && rTrack < rToGo) return 1; // going inward
4073 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4075 Int_t index=2*fCurrentEsdTrack+shieldindex;
4077 Double_t xOverX0,x0,lengthTimesMeanDensity;
4082 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4083 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4084 lengthTimesMeanDensity = xOverX0*x0;
4085 lengthTimesMeanDensity *= dir;
4086 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4089 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4090 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4093 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4094 xOverX0 = fxOverX0Shield[shieldindex];
4095 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4096 lengthTimesMeanDensity *= dir;
4097 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4100 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4101 if(fxOverX0ShieldTrks[index]<0) {
4102 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4103 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4104 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4105 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4106 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4109 xOverX0 = fxOverX0ShieldTrks[index];
4110 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4111 lengthTimesMeanDensity *= dir;
4112 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4118 //------------------------------------------------------------------------
4119 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4121 Double_t oldGlobXYZ[3],
4122 TString direction) {
4123 //-------------------------------------------------------------------
4124 // Propagate beyond layer and correct for material
4125 // (material budget in different ways according to fUseTGeo value)
4126 // Add time if going outward (PropagateTo or PropagateToTGeo)
4127 //-------------------------------------------------------------------
4129 // Define budget mode:
4130 // 0: material from AliITSRecoParam (hard coded)
4131 // 1: material from TGeo in stepsof X cm (on the fly)
4132 // X = AliITSRecoParam::GetStepSizeTGeo()
4133 // 2: material from lut
4134 // 3: material from TGeo in one step (same for all hypotheses)
4147 if(fTrackingPhase.Contains("Clusters2Tracks"))
4148 { mode=3; } else { mode=1; }
4151 if(fTrackingPhase.Contains("Clusters2Tracks"))
4152 { mode=3; } else { mode=2; }
4158 if(fTrackingPhase.Contains("Default")) mode=0;
4160 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4162 Double_t r=fgLayers[layerindex].GetR();
4163 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4165 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4167 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4169 Int_t index=6*fCurrentEsdTrack+layerindex;
4172 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4175 // back before material (no correction)
4177 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4178 if (!t->GetLocalXat(rOld,xOld)) return 0;
4179 if (!t->Propagate(xOld)) return 0;
4183 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4184 lengthTimesMeanDensity = xOverX0*x0;
4185 lengthTimesMeanDensity *= dir;
4186 // Bring the track beyond the material
4187 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4190 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4191 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4194 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4195 xOverX0 = fxOverX0Layer[layerindex];
4196 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4197 lengthTimesMeanDensity *= dir;
4198 // Bring the track beyond the material
4199 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4202 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4203 if(fxOverX0LayerTrks[index]<0) {
4204 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4205 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4206 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4207 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4208 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4209 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4212 xOverX0 = fxOverX0LayerTrks[index];
4213 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4214 lengthTimesMeanDensity *= dir;
4215 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4222 //------------------------------------------------------------------------
4223 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4224 //-----------------------------------------------------------------
4225 // Initialize LUT for storing material for each prolonged track
4226 //-----------------------------------------------------------------
4227 fxOverX0PipeTrks = new Float_t[ntracks];
4228 fxTimesRhoPipeTrks = new Float_t[ntracks];
4229 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4230 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4231 fxOverX0LayerTrks = new Float_t[ntracks*6];
4232 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4234 for(Int_t i=0; i<ntracks; i++) {
4235 fxOverX0PipeTrks[i] = -1.;
4236 fxTimesRhoPipeTrks[i] = -1.;
4238 for(Int_t j=0; j<ntracks*2; j++) {
4239 fxOverX0ShieldTrks[j] = -1.;
4240 fxTimesRhoShieldTrks[j] = -1.;
4242 for(Int_t k=0; k<ntracks*6; k++) {
4243 fxOverX0LayerTrks[k] = -1.;
4244 fxTimesRhoLayerTrks[k] = -1.;
4251 //------------------------------------------------------------------------
4252 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4253 //-----------------------------------------------------------------
4254 // Delete LUT for storing material for each prolonged track
4255 //-----------------------------------------------------------------
4256 if(fxOverX0PipeTrks) {
4257 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4259 if(fxOverX0ShieldTrks) {
4260 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4263 if(fxOverX0LayerTrks) {
4264 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4266 if(fxTimesRhoPipeTrks) {
4267 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4269 if(fxTimesRhoShieldTrks) {
4270 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4272 if(fxTimesRhoLayerTrks) {
4273 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4277 //------------------------------------------------------------------------
4278 void AliITStrackerMI::SetForceSkippingOfLayer() {
4279 //-----------------------------------------------------------------
4280 // Check if we are forced to skip layers
4281 // either we set to skip them in RecoParam
4282 // or they were off during data-taking
4283 //-----------------------------------------------------------------
4285 const AliEventInfo *eventInfo = GetEventInfo();
4287 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4288 fForceSkippingOfLayer[l] = 0;
4290 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4294 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4295 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4297 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4298 } else if(l==2 || l==3) {
4299 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4301 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4307 //------------------------------------------------------------------------
4308 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4309 Int_t ilayer,Int_t idet) const {
4310 //-----------------------------------------------------------------
4311 // This method is used to decide whether to allow a prolongation
4312 // without clusters, because we want to skip the layer.
4313 // In this case the return value is > 0:
4314 // return 1: the user requested to skip a layer
4315 // return 2: track outside z acceptance
4316 //-----------------------------------------------------------------
4318 if (ForceSkippingOfLayer(ilayer)) return 1;
4320 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4322 if (idet<0 && // out in z
4323 ilayer>innerLayCanSkip &&
4324 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4325 // check if track will cross SPD outer layer
4326 Double_t phiAtSPD2,zAtSPD2;
4327 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4328 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4330 return 2; // always allow skipping, changed on 05.11.2009
4335 //------------------------------------------------------------------------
4336 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4337 Int_t ilayer,Int_t idet,
4338 Double_t dz,Double_t dy,
4339 Bool_t noClusters) const {
4340 //-----------------------------------------------------------------
4341 // This method is used to decide whether to allow a prolongation
4342 // without clusters, because there is a dead zone in the road.
4343 // In this case the return value is > 0:
4344 // return 1: dead zone at z=0,+-7cm in SPD
4345 // return 2: all road is "bad" (dead or noisy) from the OCDB
4346 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4347 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4348 //-----------------------------------------------------------------
4350 // check dead zones at z=0,+-7cm in the SPD
4351 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4352 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4353 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4354 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4355 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4356 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4357 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4358 for (Int_t i=0; i<3; i++)
4359 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4360 AliDebug(2,Form("crack SPD %d",ilayer));
4365 // check bad zones from OCDB
4366 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4368 if (idet<0) return 0;
4370 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4373 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4374 if (ilayer==0 || ilayer==1) { // ---------- SPD
4376 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4378 detSizeFactorX *= 2.;
4379 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4382 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4383 if (detType==2) segm->SetLayer(ilayer+1);
4384 Float_t detSizeX = detSizeFactorX*segm->Dx();
4385 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4387 // check if the road overlaps with bad chips
4389 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
4390 Float_t zlocmin = zloc-dz;
4391 Float_t zlocmax = zloc+dz;
4392 Float_t xlocmin = xloc-dy;
4393 Float_t xlocmax = xloc+dy;
4394 Int_t chipsInRoad[100];
4396 // check if road goes out of detector
4397 Bool_t touchNeighbourDet=kFALSE;
4398 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4399 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4400 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4401 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4402 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));
4404 // check if this detector is bad
4406 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4407 if(!touchNeighbourDet) {
4408 return 2; // all detectors in road are bad
4410 return 3; // at least one is bad
4414 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4415 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4416 if (!nChipsInRoad) return 0;
4418 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4419 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4420 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4421 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4422 if (det.IsChipBad(chipsInRoad[iCh])) {
4430 if(!touchNeighbourDet) {
4431 AliDebug(2,"all bad in road");
4432 return 2; // all chips in road are bad
4434 return 3; // at least a bad chip in road
4439 AliDebug(2,"at least a bad in road");
4440 return 3; // at least a bad chip in road
4444 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4445 || !noClusters) return 0;
4447 // There are no clusters in road: check if there is at least
4448 // a bad SPD pixel or SDD anode or SSD strips on both sides
4450 Int_t idetInITS=idet;
4451 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4453 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4454 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4457 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4461 //------------------------------------------------------------------------
4462 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4463 const AliITStrackMI *track,
4464 Float_t &xloc,Float_t &zloc) const {
4465 //-----------------------------------------------------------------
4466 // Gives position of track in local module ref. frame
4467 //-----------------------------------------------------------------
4472 if(idet<0) return kTRUE; // track out of z acceptance of layer
4474 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4476 Int_t lad = Int_t(idet/ndet) + 1;
4478 Int_t det = idet - (lad-1)*ndet + 1;
4480 Double_t xyzGlob[3],xyzLoc[3];
4482 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4483 // take into account the misalignment: xyz at real detector plane
4484 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4486 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4488 xloc = (Float_t)xyzLoc[0];
4489 zloc = (Float_t)xyzLoc[2];
4493 //------------------------------------------------------------------------
4494 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4496 // Method to be optimized further:
4497 // Aim: decide whether a track can be used for PlaneEff evaluation
4498 // the decision is taken based on the track quality at the layer under study
4499 // no information on the clusters on this layer has to be used
4500 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4501 // the cut is done on number of sigmas from the boundaries
4503 // Input: Actual track, layer [0,5] under study
4505 // Return: kTRUE if this is a good track
4507 // it will apply a pre-selection to obtain good quality tracks.
4508 // Here also you will have the possibility to put a control on the
4509 // impact point of the track on the basic block, in order to exclude border regions
4510 // this will be done by calling a proper method of the AliITSPlaneEff class.
4512 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4513 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4515 Int_t index[AliITSgeomTGeo::kNLayers];
4517 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4519 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4520 index[k]=clusters[k];
4524 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4525 AliITSlayer &layer=fgLayers[ilayer];
4526 Double_t r=layer.GetR();
4527 AliITStrackMI tmp(*track);
4529 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4531 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
4532 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4533 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4534 if (tmp.GetClIndex(lay)>=0) ncl++;
4536 Bool_t nextout = kFALSE;
4537 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4538 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4539 Bool_t nextin = kFALSE;
4540 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4541 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4542 if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4544 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4545 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4546 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4547 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4551 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4552 Int_t idet=layer.FindDetectorIndex(phi,z);
4553 if(idet<0) { AliInfo(Form("cannot find detector"));
4556 // here check if it has good Chi Square.
4558 //propagate to the intersection with the detector plane
4559 const AliITSdetector &det=layer.GetDetector(idet);
4560 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4564 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4565 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4566 if(key>fPlaneEff->Nblock()) return kFALSE;
4567 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4568 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4570 // DEFINITION OF SEARCH ROAD FOR accepting a track
4572 //For the time being they are hard-wired, later on from AliITSRecoParam
4573 // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
4574 // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
4577 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4578 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4579 // done for RecPoints
4581 // exclude tracks at boundary between detectors
4582 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4583 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4584 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4585 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4586 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4588 if ( (locx-dx < blockXmn+boundaryWidth) ||
4589 (locx+dx > blockXmx-boundaryWidth) ||
4590 (locz-dz < blockZmn+boundaryWidth) ||
4591 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4594 //------------------------------------------------------------------------
4595 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4597 // This Method has to be optimized! For the time-being it uses the same criteria
4598 // as those used in the search of extra clusters for overlapping modules.
4600 // Method Purpose: estabilish whether a track has produced a recpoint or not
4601 // in the layer under study (For Plane efficiency)
4603 // inputs: AliITStrackMI* track (pointer to a usable track)
4605 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4606 // traversed by this very track. In details:
4607 // - if a cluster can be associated to the track then call
4608 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4609 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4612 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4613 AliITSlayer &layer=fgLayers[ilayer];
4614 Double_t r=layer.GetR();
4615 AliITStrackMI tmp(*track);
4619 if (!tmp.GetPhiZat(r,phi,z)) return;
4620 Int_t idet=layer.FindDetectorIndex(phi,z);
4622 if(idet<0) { AliInfo(Form("cannot find detector"));
4626 //propagate to the intersection with the detector plane
4627 const AliITSdetector &det=layer.GetDetector(idet);
4628 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4632 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4634 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4635 TMath::Sqrt(tmp.GetSigmaZ2() +
4636 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4637 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4638 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4639 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4640 TMath::Sqrt(tmp.GetSigmaY2() +
4641 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4642 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4643 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4645 // road in global (rphi,z) [i.e. in tracking ref. system]
4646 Double_t zmin = tmp.GetZ() - dz;
4647 Double_t zmax = tmp.GetZ() + dz;
4648 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4649 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4651 // select clusters in road
4652 layer.SelectClusters(zmin,zmax,ymin,ymax);
4654 // Define criteria for track-cluster association
4655 Double_t msz = tmp.GetSigmaZ2() +
4656 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4657 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4658 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4659 Double_t msy = tmp.GetSigmaY2() +
4660 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4661 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4662 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4663 if (tmp.GetConstrain()) {
4664 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4665 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4667 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4668 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4670 msz = 1./msz; // 1/RoadZ^2
4671 msy = 1./msy; // 1/RoadY^2
4674 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4676 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4677 //Double_t tolerance=0.2;
4678 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4679 idetc = cl->GetDetectorIndex();
4680 if(idet!=idetc) continue;
4681 //Int_t ilay = cl->GetLayer();
4683 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4684 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4686 Double_t chi2=tmp.GetPredictedChi2(cl);
4687 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4691 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4693 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4694 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4695 if(key>fPlaneEff->Nblock()) return;
4696 Bool_t found=kFALSE;
4699 while ((cl=layer.GetNextCluster(clidx))!=0) {
4700 idetc = cl->GetDetectorIndex();
4701 if(idet!=idetc) continue;
4702 // here real control to see whether the cluster can be associated to the track.
4703 // cluster not associated to track
4704 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4705 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
4706 // calculate track-clusters chi2
4707 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4708 // in particular, the error associated to the cluster
4709 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4711 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4713 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4714 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4715 // track->SetExtraModule(ilayer,idetExtra);
4717 if(!fPlaneEff->UpDatePlaneEff(found,key))
4718 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4719 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4720 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
4721 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4722 Int_t cltype[2]={-999,-999};
4726 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4727 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4730 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4731 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4732 cltype[0]=layer.GetCluster(ci)->GetNy();
4733 cltype[1]=layer.GetCluster(ci)->GetNz();
4735 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4736 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
4737 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4738 // It is computed properly by calling the method
4739 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4741 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4742 //if (tmp.PropagateTo(x,0.,0.)) {
4743 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4744 AliCluster c(*layer.GetCluster(ci));
4745 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4746 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4747 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
4748 clu[2]=TMath::Sqrt(c.GetSigmaY2());
4749 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4752 fPlaneEff->FillHistos(key,found,tr,clu,cltype);