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 "AliGeomManager.h"
38 #include "AliITSPlaneEff.h"
39 #include "AliTrackPointArray.h"
40 #include "AliESDEvent.h"
41 #include "AliESDtrack.h"
43 #include "AliITSChannelStatus.h"
44 #include "AliITSDetTypeRec.h"
45 #include "AliITSRecPoint.h"
46 #include "AliITSRecPointContainer.h"
47 #include "AliITSgeomTGeo.h"
48 #include "AliITSReconstructor.h"
49 #include "AliITSClusterParam.h"
50 #include "AliITSsegmentation.h"
51 #include "AliITSCalibration.h"
52 #include "AliITSPlaneEffSPD.h"
53 #include "AliITSPlaneEffSDD.h"
54 #include "AliITSPlaneEffSSD.h"
55 #include "AliITSV0Finder.h"
56 #include "AliITStrackerMI.h"
57 #include "AliMathBase.h"
59 ClassImp(AliITStrackerMI)
61 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
63 AliITStrackerMI::AliITStrackerMI():AliTracker(),
73 fLastLayerToTrackTo(0),
76 fTrackingPhase("Default"),
82 fxTimesRhoPipeTrks(0),
83 fxOverX0ShieldTrks(0),
84 fxTimesRhoShieldTrks(0),
86 fxTimesRhoLayerTrks(0),
93 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
94 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
95 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
97 //------------------------------------------------------------------------
98 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
99 fI(AliITSgeomTGeo::GetNLayers()),
108 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
111 fTrackingPhase("Default"),
117 fxTimesRhoPipeTrks(0),
118 fxOverX0ShieldTrks(0),
119 fxTimesRhoShieldTrks(0),
120 fxOverX0LayerTrks(0),
121 fxTimesRhoLayerTrks(0),
123 fITSChannelStatus(0),
126 //--------------------------------------------------------------------
127 //This is the AliITStrackerMI constructor
128 //--------------------------------------------------------------------
130 AliWarning("\"geom\" is actually a dummy argument !");
136 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
137 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
138 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
140 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
141 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
142 Double_t poff=TMath::ATan2(y,x);
145 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
146 Double_t r=TMath::Sqrt(x*x + y*y);
148 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
149 r += TMath::Sqrt(x*x + y*y);
150 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
151 r += TMath::Sqrt(x*x + y*y);
152 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
153 r += TMath::Sqrt(x*x + y*y);
156 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
158 for (Int_t j=1; j<nlad+1; j++) {
159 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
160 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
161 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
163 Double_t txyz[3]={0.};
164 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
165 m.LocalToMaster(txyz,xyz);
166 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
167 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
169 if (phi<0) phi+=TMath::TwoPi();
170 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
172 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
173 new(&det) AliITSdetector(r,phi);
174 // compute the real radius (with misalignment)
175 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
177 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
178 mmisal.LocalToMaster(txyz,xyz);
179 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
180 det.SetRmisal(rmisal);
182 } // end loop on detectors
183 } // end loop on ladders
184 fForceSkippingOfLayer[i] = 0;
185 } // end loop on layers
188 fI=AliITSgeomTGeo::GetNLayers();
191 fConstraint[0]=1; fConstraint[1]=0;
193 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
194 AliITSReconstructor::GetRecoParam()->GetYVdef(),
195 AliITSReconstructor::GetRecoParam()->GetZVdef()};
196 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
197 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
198 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
199 SetVertex(xyzVtx,ersVtx);
201 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
202 for (Int_t i=0;i<100000;i++){
203 fBestTrackIndex[i]=0;
206 // store positions of centre of SPD modules (in z)
207 // The convetion is that fSPDdetzcentre is ordered from -z to +z
209 if (tr[2]<0) { // old geom (up to v5asymmPPR)
210 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
211 fSPDdetzcentre[0] = tr[2];
212 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
213 fSPDdetzcentre[1] = tr[2];
214 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
215 fSPDdetzcentre[2] = tr[2];
216 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
217 fSPDdetzcentre[3] = tr[2];
218 } else { // new geom (from v11Hybrid)
219 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
220 fSPDdetzcentre[0] = tr[2];
221 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
222 fSPDdetzcentre[1] = tr[2];
223 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
224 fSPDdetzcentre[2] = tr[2];
225 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
226 fSPDdetzcentre[3] = tr[2];
229 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
230 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
231 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
235 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
236 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
238 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
240 // only for plane efficiency evaluation
241 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
242 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
243 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
244 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
245 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
246 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
247 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
248 else fPlaneEff = new AliITSPlaneEffSSD();
249 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
250 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
251 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
254 //------------------------------------------------------------------------
255 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
257 fBestTrack(tracker.fBestTrack),
258 fTrackToFollow(tracker.fTrackToFollow),
259 fTrackHypothesys(tracker.fTrackHypothesys),
260 fBestHypothesys(tracker.fBestHypothesys),
261 fOriginal(tracker.fOriginal),
262 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
263 fPass(tracker.fPass),
264 fAfterV0(tracker.fAfterV0),
265 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
266 fCoefficients(tracker.fCoefficients),
268 fTrackingPhase(tracker.fTrackingPhase),
269 fUseTGeo(tracker.fUseTGeo),
270 fNtracks(tracker.fNtracks),
271 fxOverX0Pipe(tracker.fxOverX0Pipe),
272 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
274 fxTimesRhoPipeTrks(0),
275 fxOverX0ShieldTrks(0),
276 fxTimesRhoShieldTrks(0),
277 fxOverX0LayerTrks(0),
278 fxTimesRhoLayerTrks(0),
279 fDebugStreamer(tracker.fDebugStreamer),
280 fITSChannelStatus(tracker.fITSChannelStatus),
281 fkDetTypeRec(tracker.fkDetTypeRec),
282 fPlaneEff(tracker.fPlaneEff) {
286 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
289 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
290 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
293 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
294 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
297 //------------------------------------------------------------------------
298 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
299 //Assignment operator
300 this->~AliITStrackerMI();
301 new(this) AliITStrackerMI(tracker);
304 //------------------------------------------------------------------------
305 AliITStrackerMI::~AliITStrackerMI()
310 if (fCoefficients) delete [] fCoefficients;
311 DeleteTrksMaterialLUT();
312 if (fDebugStreamer) {
313 //fDebugStreamer->Close();
314 delete fDebugStreamer;
316 if(fITSChannelStatus) delete fITSChannelStatus;
317 if(fPlaneEff) delete fPlaneEff;
319 //------------------------------------------------------------------------
320 void AliITStrackerMI::ReadBadFromDetTypeRec() {
321 //--------------------------------------------------------------------
322 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
324 //--------------------------------------------------------------------
326 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
328 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
330 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
333 if(fITSChannelStatus) delete fITSChannelStatus;
334 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
336 // ITS detectors and chips
337 Int_t i=0,j=0,k=0,ndet=0;
338 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
339 Int_t nBadDetsPerLayer=0;
340 ndet=AliITSgeomTGeo::GetNDetectors(i);
341 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
342 for (k=1; k<ndet+1; k++) {
343 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
344 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
345 if(det.IsBad()) {nBadDetsPerLayer++;}
346 } // end loop on detectors
347 } // end loop on ladders
348 Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
349 } // end loop on layers
353 //------------------------------------------------------------------------
354 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
355 //--------------------------------------------------------------------
356 //This function loads ITS clusters
357 //--------------------------------------------------------------------
359 TClonesArray *clusters = NULL;
360 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
361 clusters=rpcont->FetchClusters(0,cTree);
362 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
363 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
366 Int_t i=0,j=0,ndet=0;
368 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
369 ndet=fgLayers[i].GetNdetectors();
370 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
371 for (; j<jmax; j++) {
372 // if (!cTree->GetEvent(j)) continue;
373 clusters = rpcont->UncheckedGetClusters(j);
374 if(!clusters)continue;
375 Int_t ncl=clusters->GetEntriesFast();
376 SignDeltas(clusters,GetZ());
379 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
380 detector=c->GetDetectorIndex();
382 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
384 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
386 AliWarning(Form("Too many clusters on layer %d!",i));
391 // add dead zone "virtual" cluster in SPD, if there is a cluster within
392 // zwindow cm from the dead zone
393 // This method assumes that fSPDdetzcentre is ordered from -z to +z
394 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
395 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
396 Int_t lab[4] = {0,0,0,detector};
397 Int_t info[3] = {0,0,i};
398 Float_t q = 0.; // this identifies virtual clusters
399 Float_t hit[5] = {xdead,
401 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
402 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
404 Bool_t local = kTRUE;
405 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
406 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
407 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
408 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
409 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
410 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
411 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
412 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
413 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
414 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
415 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
416 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
417 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
418 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
419 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
420 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
421 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
422 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
423 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
425 } // "virtual" clusters in SPD
429 fgLayers[i].ResetRoad(); //road defined by the cluster density
430 fgLayers[i].SortClusters();
433 // check whether we have to skip some layers
434 SetForceSkippingOfLayer();
438 //------------------------------------------------------------------------
439 void AliITStrackerMI::UnloadClusters() {
440 //--------------------------------------------------------------------
441 //This function unloads ITS clusters
442 //--------------------------------------------------------------------
443 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
445 //------------------------------------------------------------------------
446 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
447 //--------------------------------------------------------------------
448 // Publishes all pointers to clusters known to the tracker into the
449 // passed object array.
450 // The ownership is not transfered - the caller is not expected to delete
452 //--------------------------------------------------------------------
454 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
455 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
456 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
463 //------------------------------------------------------------------------
464 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
465 //--------------------------------------------------------------------
466 // Correction for the material between the TPC and the ITS
467 //--------------------------------------------------------------------
468 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
469 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
470 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
471 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
472 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
473 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
474 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
475 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
477 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
483 //------------------------------------------------------------------------
484 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
485 //--------------------------------------------------------------------
486 // This functions reconstructs ITS tracks
487 // The clusters must be already loaded !
488 //--------------------------------------------------------------------
490 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
492 fTrackingPhase="Clusters2Tracks";
494 TObjArray itsTracks(15000);
496 fEsd = event; // store pointer to the esd
498 // temporary (for cosmics)
499 if(event->GetVertex()) {
500 TString title = event->GetVertex()->GetTitle();
501 if(title.Contains("cosmics")) {
502 Double_t xyz[3]={GetX(),GetY(),GetZ()};
503 Double_t exyz[3]={0.1,0.1,0.1};
509 {/* Read ESD tracks */
510 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
511 Int_t nentr=event->GetNumberOfTracks();
513 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
515 AliESDtrack *esd=event->GetTrack(nentr);
516 // ---- for debugging:
517 //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); }
519 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
520 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
521 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
522 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
525 t=new AliITStrackMI(*esd);
526 } catch (const Char_t *msg) {
527 //Warning("Clusters2Tracks",msg);
531 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
532 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
535 // look at the ESD mass hypothesys !
536 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
538 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
540 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
541 //track - can be V0 according to TPC
543 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
547 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
551 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
555 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
560 t->SetReconstructed(kFALSE);
561 itsTracks.AddLast(t);
562 fOriginal.AddLast(t);
564 } /* End Read ESD tracks */
568 Int_t nentr=itsTracks.GetEntriesFast();
569 fTrackHypothesys.Expand(nentr);
570 fBestHypothesys.Expand(nentr);
571 MakeCoefficients(nentr);
572 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
574 // THE TWO TRACKING PASSES
575 for (fPass=0; fPass<2; fPass++) {
576 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
577 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
578 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
579 if (t==0) continue; //this track has been already tracked
580 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
581 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
582 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
583 if (fConstraint[fPass]) {
584 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
585 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
588 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
589 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
591 ResetTrackToFollow(*t);
594 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
597 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
599 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
600 if (!besttrack) continue;
601 besttrack->SetLabel(tpcLabel);
602 // besttrack->CookdEdx();
604 besttrack->SetFakeRatio(1.);
605 CookLabel(besttrack,0.); //For comparison only
606 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
608 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
610 t->SetReconstructed(kTRUE);
612 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
614 GetBestHypothesysMIP(itsTracks);
615 } // end loop on the two tracking passes
617 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
618 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
623 Int_t entries = fTrackHypothesys.GetEntriesFast();
624 for (Int_t ientry=0; ientry<entries; ientry++) {
625 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
626 if (array) array->Delete();
627 delete fTrackHypothesys.RemoveAt(ientry);
630 fTrackHypothesys.Delete();
631 fBestHypothesys.Delete();
633 delete [] fCoefficients;
635 DeleteTrksMaterialLUT();
637 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
639 fTrackingPhase="Default";
643 //------------------------------------------------------------------------
644 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
645 //--------------------------------------------------------------------
646 // This functions propagates reconstructed ITS tracks back
647 // The clusters must be loaded !
648 //--------------------------------------------------------------------
649 fTrackingPhase="PropagateBack";
650 Int_t nentr=event->GetNumberOfTracks();
651 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
654 for (Int_t i=0; i<nentr; i++) {
655 AliESDtrack *esd=event->GetTrack(i);
657 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
658 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
662 t=new AliITStrackMI(*esd);
663 } catch (const Char_t *msg) {
664 //Warning("PropagateBack",msg);
668 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
670 ResetTrackToFollow(*t);
673 // propagate to vertex [SR, GSI 17.02.2003]
674 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
675 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
676 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
677 fTrackToFollow.StartTimeIntegral();
678 // from vertex to outside pipe
679 CorrectForPipeMaterial(&fTrackToFollow,"outward");
681 // Start time integral and add distance from current position to vertex
682 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
683 fTrackToFollow.GetXYZ(xyzTrk);
685 for (Int_t icoord=0; icoord<3; icoord++) {
686 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
689 fTrackToFollow.StartTimeIntegral();
690 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
692 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
693 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
694 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
698 fTrackToFollow.SetLabel(t->GetLabel());
699 //fTrackToFollow.CookdEdx();
700 CookLabel(&fTrackToFollow,0.); //For comparison only
701 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
702 //UseClusters(&fTrackToFollow);
708 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
710 fTrackingPhase="Default";
714 //------------------------------------------------------------------------
715 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
716 //--------------------------------------------------------------------
717 // This functions refits ITS tracks using the
718 // "inward propagated" TPC tracks
719 // The clusters must be loaded !
720 //--------------------------------------------------------------------
721 fTrackingPhase="RefitInward";
723 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
725 Int_t nentr=event->GetNumberOfTracks();
726 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
729 for (Int_t i=0; i<nentr; i++) {
730 AliESDtrack *esd=event->GetTrack(i);
732 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
733 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
734 if (esd->GetStatus()&AliESDtrack::kTPCout)
735 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
739 t=new AliITStrackMI(*esd);
740 } catch (const Char_t *msg) {
741 //Warning("RefitInward",msg);
745 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
746 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
751 ResetTrackToFollow(*t);
752 fTrackToFollow.ResetClusters();
754 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
755 fTrackToFollow.ResetCovariance(10.);
758 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
759 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
761 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
762 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
763 AliDebug(2," refit OK");
764 fTrackToFollow.SetLabel(t->GetLabel());
765 // fTrackToFollow.CookdEdx();
766 CookdEdx(&fTrackToFollow);
768 CookLabel(&fTrackToFollow,0.0); //For comparison only
771 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
772 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
773 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
774 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
775 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
776 Double_t r[3]={0.,0.,0.};
778 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
785 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
787 fTrackingPhase="Default";
791 //------------------------------------------------------------------------
792 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
793 //--------------------------------------------------------------------
794 // Return pointer to a given cluster
795 //--------------------------------------------------------------------
796 Int_t l=(index & 0xf0000000) >> 28;
797 Int_t c=(index & 0x0fffffff) >> 00;
798 return fgLayers[l].GetCluster(c);
800 //------------------------------------------------------------------------
801 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
802 //--------------------------------------------------------------------
803 // Get track space point with index i
804 //--------------------------------------------------------------------
806 Int_t l=(index & 0xf0000000) >> 28;
807 Int_t c=(index & 0x0fffffff) >> 00;
808 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
809 Int_t idet = cl->GetDetectorIndex();
813 cl->GetGlobalXYZ(xyz);
814 cl->GetGlobalCov(cov);
816 p.SetCharge(cl->GetQ());
817 p.SetDriftTime(cl->GetDriftTime());
818 p.SetChargeRatio(cl->GetChargeRatio());
819 p.SetClusterType(cl->GetClusterType());
820 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
823 iLayer = AliGeomManager::kSPD1;
826 iLayer = AliGeomManager::kSPD2;
829 iLayer = AliGeomManager::kSDD1;
832 iLayer = AliGeomManager::kSDD2;
835 iLayer = AliGeomManager::kSSD1;
838 iLayer = AliGeomManager::kSSD2;
841 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
844 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
845 p.SetVolumeID((UShort_t)volid);
848 //------------------------------------------------------------------------
849 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
850 AliTrackPoint& p, const AliESDtrack *t) {
851 //--------------------------------------------------------------------
852 // Get track space point with index i
853 // (assign error estimated during the tracking)
854 //--------------------------------------------------------------------
856 Int_t l=(index & 0xf0000000) >> 28;
857 Int_t c=(index & 0x0fffffff) >> 00;
858 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
859 Int_t idet = cl->GetDetectorIndex();
861 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
863 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
865 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
866 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
867 Double_t alpha = t->GetAlpha();
868 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
869 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
870 phi += alpha-det.GetPhi();
871 Float_t tgphi = TMath::Tan(phi);
873 Float_t tgl = t->GetTgl(); // tgl about const along track
874 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
876 Float_t errtrky,errtrkz,covyz;
877 Bool_t addMisalErr=kFALSE;
878 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
882 cl->GetGlobalXYZ(xyz);
883 // cl->GetGlobalCov(cov);
884 Float_t pos[3] = {0.,0.,0.};
885 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
886 tmpcl.GetGlobalCov(cov);
889 p.SetCharge(cl->GetQ());
890 p.SetDriftTime(cl->GetDriftTime());
891 p.SetChargeRatio(cl->GetChargeRatio());
892 p.SetClusterType(cl->GetClusterType());
894 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
897 iLayer = AliGeomManager::kSPD1;
900 iLayer = AliGeomManager::kSPD2;
903 iLayer = AliGeomManager::kSDD1;
906 iLayer = AliGeomManager::kSDD2;
909 iLayer = AliGeomManager::kSSD1;
912 iLayer = AliGeomManager::kSSD2;
915 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
918 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
920 p.SetVolumeID((UShort_t)volid);
923 //------------------------------------------------------------------------
924 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
926 //--------------------------------------------------------------------
927 // Follow prolongation tree
928 //--------------------------------------------------------------------
930 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
931 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
934 AliESDtrack * esd = otrack->GetESDtrack();
935 if (esd->GetV0Index(0)>0) {
936 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
937 // mapping of ESD track is different as ITS track in Containers
938 // Need something more stable
939 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
940 for (Int_t i=0;i<3;i++){
941 Int_t index = esd->GetV0Index(i);
943 AliESDv0 * vertex = fEsd->GetV0(index);
944 if (vertex->GetStatus()<0) continue; // rejected V0
946 if (esd->GetSign()>0) {
947 vertex->SetIndex(0,esdindex);
949 vertex->SetIndex(1,esdindex);
953 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
955 bestarray = new TObjArray(5);
956 fBestHypothesys.AddAt(bestarray,esdindex);
960 //setup tree of the prolongations
962 static AliITStrackMI tracks[7][100];
963 AliITStrackMI *currenttrack;
964 static AliITStrackMI currenttrack1;
965 static AliITStrackMI currenttrack2;
966 static AliITStrackMI backuptrack;
968 Int_t nindexes[7][100];
969 Float_t normalizedchi2[100];
970 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
971 otrack->SetNSkipped(0);
972 new (&(tracks[6][0])) AliITStrackMI(*otrack);
974 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
975 Int_t modstatus = 1; // found
979 // follow prolongations
980 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
981 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
984 AliITSlayer &layer=fgLayers[ilayer];
985 Double_t r = layer.GetR();
991 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
993 if (ntracks[ilayer]>=100) break;
994 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
995 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
996 if (ntracks[ilayer]>15+ilayer){
997 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
998 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1001 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1003 // material between SSD and SDD, SDD and SPD
1005 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1007 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1011 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1012 Int_t idet=layer.FindDetectorIndex(phi,z);
1014 Double_t trackGlobXYZ1[3];
1015 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1017 // Get the budget to the primary vertex for the current track being prolonged
1018 Double_t budgetToPrimVertex = GetEffectiveThickness();
1020 // check if we allow a prolongation without point
1021 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1023 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1024 // propagate to the layer radius
1025 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1026 if(!vtrack->Propagate(xToGo)) continue;
1027 // apply correction for material of the current layer
1028 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1029 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1030 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1031 vtrack->SetClIndex(ilayer,-1);
1032 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1033 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1034 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1036 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1041 // track outside layer acceptance in z
1042 if (idet<0) continue;
1044 //propagate to the intersection with the detector plane
1045 const AliITSdetector &det=layer.GetDetector(idet);
1046 new(¤ttrack2) AliITStrackMI(currenttrack1);
1047 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1048 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1049 currenttrack1.SetDetectorIndex(idet);
1050 currenttrack2.SetDetectorIndex(idet);
1051 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1054 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1056 // road in global (rphi,z) [i.e. in tracking ref. system]
1057 Double_t zmin,zmax,ymin,ymax;
1058 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1060 // select clusters in road
1061 layer.SelectClusters(zmin,zmax,ymin,ymax);
1062 //********************
1064 // Define criteria for track-cluster association
1065 Double_t msz = currenttrack1.GetSigmaZ2() +
1066 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1067 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1068 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1069 Double_t msy = currenttrack1.GetSigmaY2() +
1070 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1071 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1072 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1074 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1075 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1077 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1078 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1080 msz = 1./msz; // 1/RoadZ^2
1081 msy = 1./msy; // 1/RoadY^2
1085 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1087 const AliITSRecPoint *cl=0;
1089 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1090 Bool_t deadzoneSPD=kFALSE;
1091 currenttrack = ¤ttrack1;
1093 // check if the road contains a dead zone
1094 Bool_t noClusters = kFALSE;
1095 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1096 if (noClusters) AliDebug(2,"no clusters in road");
1097 Double_t dz=0.5*(zmax-zmin);
1098 Double_t dy=0.5*(ymax-ymin);
1099 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1100 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1101 // create a prolongation without clusters (check also if there are no clusters in the road)
1104 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1105 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1106 updatetrack->SetClIndex(ilayer,-1);
1108 modstatus = 5; // no cls in road
1109 } else if (dead==1) {
1110 modstatus = 7; // holes in z in SPD
1111 } else if (dead==2 || dead==3 || dead==4) {
1112 modstatus = 2; // dead from OCDB
1114 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1115 // apply correction for material of the current layer
1116 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1117 if (constrain) { // apply vertex constrain
1118 updatetrack->SetConstrain(constrain);
1119 Bool_t isPrim = kTRUE;
1120 if (ilayer<4) { // check that it's close to the vertex
1121 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1122 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1123 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1124 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1125 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1127 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1129 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1131 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1132 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1134 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1135 updatetrack->SetDeadZoneProbability(ilayer,1.);
1136 } else if (dead==4) { // at least a single dead channel from OCDB
1137 updatetrack->SetDeadZoneProbability(ilayer,0.);
1144 // loop over clusters in the road
1145 while ((cl=layer.GetNextCluster(clidx))!=0) {
1146 if (ntracks[ilayer]>95) break; //space for skipped clusters
1147 Bool_t changedet =kFALSE;
1148 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1149 Int_t idetc=cl->GetDetectorIndex();
1151 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1152 // take into account misalignment (bring track to real detector plane)
1153 Double_t xTrOrig = currenttrack->GetX();
1154 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1155 // a first cut on track-cluster distance
1156 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1157 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1158 { // cluster not associated to track
1159 AliDebug(2,"not associated");
1162 // bring track back to ideal detector plane
1163 if (!currenttrack->Propagate(xTrOrig)) continue;
1164 } else { // have to move track to cluster's detector
1165 const AliITSdetector &detc=layer.GetDetector(idetc);
1166 // a first cut on track-cluster distance
1168 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1169 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1170 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1171 continue; // cluster not associated to track
1173 new (&backuptrack) AliITStrackMI(currenttrack2);
1175 currenttrack =¤ttrack2;
1176 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1177 new (currenttrack) AliITStrackMI(backuptrack);
1181 currenttrack->SetDetectorIndex(idetc);
1182 // Get again the budget to the primary vertex
1183 // for the current track being prolonged, if had to change detector
1184 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1187 // calculate track-clusters chi2
1188 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1190 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1191 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1192 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1193 if (ntracks[ilayer]>=100) continue;
1194 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1195 updatetrack->SetClIndex(ilayer,-1);
1196 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1198 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1199 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1200 AliDebug(2,"update failed");
1203 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1204 modstatus = 1; // found
1205 } else { // virtual cluster in dead zone
1206 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1207 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1208 modstatus = 7; // holes in z in SPD
1212 Float_t xlocnewdet,zlocnewdet;
1213 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1214 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1217 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1219 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1221 // apply correction for material of the current layer
1222 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1224 if (constrain) { // apply vertex constrain
1225 updatetrack->SetConstrain(constrain);
1226 Bool_t isPrim = kTRUE;
1227 if (ilayer<4) { // check that it's close to the vertex
1228 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1229 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1230 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1231 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1232 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1234 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1235 } //apply vertex constrain
1237 } // create new hypothesis
1239 AliDebug(2,"chi2 too large");
1242 } // loop over possible prolongations
1244 // allow one prolongation without clusters
1245 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1246 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1247 // apply correction for material of the current layer
1248 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1249 vtrack->SetClIndex(ilayer,-1);
1250 modstatus = 3; // skipped
1251 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1252 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1253 vtrack->IncrementNSkipped();
1258 } // loop over tracks in layer ilayer+1
1260 //loop over track candidates for the current layer
1266 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1267 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1268 if (normalizedchi2[itrack] <
1269 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1273 if (constrain) { // constrain
1274 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1276 } else { // !constrain
1277 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1282 // sort tracks by increasing normalized chi2
1283 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1284 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1285 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1286 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1287 } // end loop over layers
1291 // Now select tracks to be kept
1293 Int_t max = constrain ? 20 : 5;
1295 // tracks that reach layer 0 (SPD inner)
1296 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1297 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1298 if (track.GetNumberOfClusters()<2) continue;
1299 if (!constrain && track.GetNormChi2(0) >
1300 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1303 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1306 // tracks that reach layer 1 (SPD outer)
1307 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1308 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1309 if (track.GetNumberOfClusters()<4) continue;
1310 if (!constrain && track.GetNormChi2(1) >
1311 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1312 if (constrain) track.IncrementNSkipped();
1314 track.SetD(0,track.GetD(GetX(),GetY()));
1315 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1316 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1317 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1320 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1323 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1325 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1326 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1327 if (track.GetNumberOfClusters()<3) continue;
1328 if (!constrain && track.GetNormChi2(2) >
1329 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1330 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1332 track.SetD(0,track.GetD(GetX(),GetY()));
1333 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1334 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1335 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1338 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1344 // register best track of each layer - important for V0 finder
1346 for (Int_t ilayer=0;ilayer<5;ilayer++){
1347 if (ntracks[ilayer]==0) continue;
1348 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1349 if (track.GetNumberOfClusters()<1) continue;
1350 CookLabel(&track,0);
1351 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1355 // update TPC V0 information
1357 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1358 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1359 for (Int_t i=0;i<3;i++){
1360 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1361 if (index==0) break;
1362 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1363 if (vertex->GetStatus()<0) continue; // rejected V0
1365 if (otrack->GetSign()>0) {
1366 vertex->SetIndex(0,esdindex);
1369 vertex->SetIndex(1,esdindex);
1371 //find nearest layer with track info
1372 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1373 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1374 Int_t nearest = nearestold;
1375 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1376 if (ntracks[nearest]==0){
1381 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1382 if (nearestold<5&&nearest<5){
1383 Bool_t accept = track.GetNormChi2(nearest)<10;
1385 if (track.GetSign()>0) {
1386 vertex->SetParamP(track);
1387 vertex->Update(fprimvertex);
1388 //vertex->SetIndex(0,track.fESDtrack->GetID());
1389 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1391 vertex->SetParamN(track);
1392 vertex->Update(fprimvertex);
1393 //vertex->SetIndex(1,track.fESDtrack->GetID());
1394 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1396 vertex->SetStatus(vertex->GetStatus()+1);
1398 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1405 //------------------------------------------------------------------------
1406 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1408 //--------------------------------------------------------------------
1411 return fgLayers[layer];
1413 //------------------------------------------------------------------------
1414 AliITStrackerMI::AliITSlayer::AliITSlayer():
1444 //--------------------------------------------------------------------
1445 //default AliITSlayer constructor
1446 //--------------------------------------------------------------------
1447 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1448 fClusterWeight[i]=0;
1449 fClusterTracks[0][i]=-1;
1450 fClusterTracks[1][i]=-1;
1451 fClusterTracks[2][i]=-1;
1452 fClusterTracks[3][i]=-1;
1455 //------------------------------------------------------------------------
1456 AliITStrackerMI::AliITSlayer::
1457 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1486 //--------------------------------------------------------------------
1487 //main AliITSlayer constructor
1488 //--------------------------------------------------------------------
1489 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1490 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1492 //------------------------------------------------------------------------
1493 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1495 fPhiOffset(layer.fPhiOffset),
1496 fNladders(layer.fNladders),
1497 fZOffset(layer.fZOffset),
1498 fNdetectors(layer.fNdetectors),
1499 fDetectors(layer.fDetectors),
1504 fClustersCs(layer.fClustersCs),
1505 fClusterIndexCs(layer.fClusterIndexCs),
1509 fCurrentSlice(layer.fCurrentSlice),
1517 fAccepted(layer.fAccepted),
1519 fMaxSigmaClY(layer.fMaxSigmaClY),
1520 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1521 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1525 //------------------------------------------------------------------------
1526 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1527 //--------------------------------------------------------------------
1528 // AliITSlayer destructor
1529 //--------------------------------------------------------------------
1530 delete [] fDetectors;
1531 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1532 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1533 fClusterWeight[i]=0;
1534 fClusterTracks[0][i]=-1;
1535 fClusterTracks[1][i]=-1;
1536 fClusterTracks[2][i]=-1;
1537 fClusterTracks[3][i]=-1;
1540 //------------------------------------------------------------------------
1541 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1542 //--------------------------------------------------------------------
1543 // This function removes loaded clusters
1544 //--------------------------------------------------------------------
1545 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1546 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1547 fClusterWeight[i]=0;
1548 fClusterTracks[0][i]=-1;
1549 fClusterTracks[1][i]=-1;
1550 fClusterTracks[2][i]=-1;
1551 fClusterTracks[3][i]=-1;
1557 //------------------------------------------------------------------------
1558 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1559 //--------------------------------------------------------------------
1560 // This function reset weights of the clusters
1561 //--------------------------------------------------------------------
1562 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1563 fClusterWeight[i]=0;
1564 fClusterTracks[0][i]=-1;
1565 fClusterTracks[1][i]=-1;
1566 fClusterTracks[2][i]=-1;
1567 fClusterTracks[3][i]=-1;
1569 for (Int_t i=0; i<fN;i++) {
1570 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1571 if (cl&&cl->IsUsed()) cl->Use();
1575 //------------------------------------------------------------------------
1576 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1577 //--------------------------------------------------------------------
1578 // This function calculates the road defined by the cluster density
1579 //--------------------------------------------------------------------
1581 for (Int_t i=0; i<fN; i++) {
1582 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1584 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1586 //------------------------------------------------------------------------
1587 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1588 //--------------------------------------------------------------------
1589 //This function adds a cluster to this layer
1590 //--------------------------------------------------------------------
1591 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1597 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1599 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1600 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1601 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1602 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1603 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1604 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1607 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1608 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1609 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1610 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1614 //------------------------------------------------------------------------
1615 void AliITStrackerMI::AliITSlayer::SortClusters()
1620 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1621 Float_t *z = new Float_t[fN];
1622 Int_t * index = new Int_t[fN];
1624 fMaxSigmaClY=0.; //AD
1625 fMaxSigmaClZ=0.; //AD
1627 for (Int_t i=0;i<fN;i++){
1628 z[i] = fClusters[i]->GetZ();
1629 // save largest errors in y and z for this layer
1630 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1631 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1633 TMath::Sort(fN,z,index,kFALSE);
1634 for (Int_t i=0;i<fN;i++){
1635 clusters[i] = fClusters[index[i]];
1638 for (Int_t i=0;i<fN;i++){
1639 fClusters[i] = clusters[i];
1640 fZ[i] = fClusters[i]->GetZ();
1641 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1642 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1643 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1653 for (Int_t i=0;i<fN;i++){
1654 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1655 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1656 fClusterIndex[i] = i;
1660 fDy5 = (fYB[1]-fYB[0])/5.;
1661 fDy10 = (fYB[1]-fYB[0])/10.;
1662 fDy20 = (fYB[1]-fYB[0])/20.;
1663 for (Int_t i=0;i<6;i++) fN5[i] =0;
1664 for (Int_t i=0;i<11;i++) fN10[i]=0;
1665 for (Int_t i=0;i<21;i++) fN20[i]=0;
1667 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;}
1668 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;}
1669 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;}
1672 for (Int_t i=0;i<fN;i++)
1673 for (Int_t irot=-1;irot<=1;irot++){
1674 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1676 for (Int_t slice=0; slice<6;slice++){
1677 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1678 fClusters5[slice][fN5[slice]] = fClusters[i];
1679 fY5[slice][fN5[slice]] = curY;
1680 fZ5[slice][fN5[slice]] = fZ[i];
1681 fClusterIndex5[slice][fN5[slice]]=i;
1686 for (Int_t slice=0; slice<11;slice++){
1687 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1688 fClusters10[slice][fN10[slice]] = fClusters[i];
1689 fY10[slice][fN10[slice]] = curY;
1690 fZ10[slice][fN10[slice]] = fZ[i];
1691 fClusterIndex10[slice][fN10[slice]]=i;
1696 for (Int_t slice=0; slice<21;slice++){
1697 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1698 fClusters20[slice][fN20[slice]] = fClusters[i];
1699 fY20[slice][fN20[slice]] = curY;
1700 fZ20[slice][fN20[slice]] = fZ[i];
1701 fClusterIndex20[slice][fN20[slice]]=i;
1708 // consistency check
1710 for (Int_t i=0;i<fN-1;i++){
1716 for (Int_t slice=0;slice<21;slice++)
1717 for (Int_t i=0;i<fN20[slice]-1;i++){
1718 if (fZ20[slice][i]>fZ20[slice][i+1]){
1725 //------------------------------------------------------------------------
1726 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1727 //--------------------------------------------------------------------
1728 // This function returns the index of the nearest cluster
1729 //--------------------------------------------------------------------
1732 if (fCurrentSlice<0) {
1741 if (ncl==0) return 0;
1742 Int_t b=0, e=ncl-1, m=(b+e)/2;
1743 for (; b<e; m=(b+e)/2) {
1744 // if (z > fClusters[m]->GetZ()) b=m+1;
1745 if (z > zcl[m]) b=m+1;
1750 //------------------------------------------------------------------------
1751 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 {
1752 //--------------------------------------------------------------------
1753 // This function computes the rectangular road for this track
1754 //--------------------------------------------------------------------
1757 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1758 // take into account the misalignment: propagate track to misaligned detector plane
1759 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1761 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1762 TMath::Sqrt(track->GetSigmaZ2() +
1763 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1764 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1765 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1766 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1767 TMath::Sqrt(track->GetSigmaY2() +
1768 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1769 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1770 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1772 // track at boundary between detectors, enlarge road
1773 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1774 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1775 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1776 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1777 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1778 Float_t tgl = TMath::Abs(track->GetTgl());
1779 if (tgl > 1.) tgl=1.;
1780 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1781 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1782 Float_t snp = TMath::Abs(track->GetSnp());
1783 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1784 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1787 // add to the road a term (up to 2-3 mm) to deal with misalignments
1788 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1789 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1791 Double_t r = fgLayers[ilayer].GetR();
1792 zmin = track->GetZ() - dz;
1793 zmax = track->GetZ() + dz;
1794 ymin = track->GetY() + r*det.GetPhi() - dy;
1795 ymax = track->GetY() + r*det.GetPhi() + dy;
1797 // bring track back to idead detector plane
1798 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1802 //------------------------------------------------------------------------
1803 void AliITStrackerMI::AliITSlayer::
1804 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1805 //--------------------------------------------------------------------
1806 // This function sets the "window"
1807 //--------------------------------------------------------------------
1809 Double_t circle=2*TMath::Pi()*fR;
1815 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1816 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1817 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1818 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1819 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1821 Float_t ymiddle = (fYmin+fYmax)*0.5;
1822 if (ymiddle<fYB[0]) {
1823 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1824 } else if (ymiddle>fYB[1]) {
1825 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1831 fClustersCs = fClusters;
1832 fClusterIndexCs = fClusterIndex;
1838 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1839 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1840 if (slice<0) slice=0;
1841 if (slice>20) slice=20;
1842 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1844 fCurrentSlice=slice;
1845 fClustersCs = fClusters20[fCurrentSlice];
1846 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1847 fYcs = fY20[fCurrentSlice];
1848 fZcs = fZ20[fCurrentSlice];
1849 fNcs = fN20[fCurrentSlice];
1854 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1855 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1856 if (slice<0) slice=0;
1857 if (slice>10) slice=10;
1858 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1860 fCurrentSlice=slice;
1861 fClustersCs = fClusters10[fCurrentSlice];
1862 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1863 fYcs = fY10[fCurrentSlice];
1864 fZcs = fZ10[fCurrentSlice];
1865 fNcs = fN10[fCurrentSlice];
1870 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1871 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1872 if (slice<0) slice=0;
1873 if (slice>5) slice=5;
1874 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1876 fCurrentSlice=slice;
1877 fClustersCs = fClusters5[fCurrentSlice];
1878 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1879 fYcs = fY5[fCurrentSlice];
1880 fZcs = fZ5[fCurrentSlice];
1881 fNcs = fN5[fCurrentSlice];
1885 fI = FindClusterIndex(fZmin);
1886 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
1892 //------------------------------------------------------------------------
1893 Int_t AliITStrackerMI::AliITSlayer::
1894 FindDetectorIndex(Double_t phi, Double_t z) const {
1895 //--------------------------------------------------------------------
1896 //This function finds the detector crossed by the track
1897 //--------------------------------------------------------------------
1899 if (fZOffset<0) // old geometry
1900 dphi = -(phi-fPhiOffset);
1901 else // new geometry
1902 dphi = phi-fPhiOffset;
1905 if (dphi < 0) dphi += 2*TMath::Pi();
1906 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1907 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1908 if (np>=fNladders) np-=fNladders;
1909 if (np<0) np+=fNladders;
1912 Double_t dz=fZOffset-z;
1913 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1914 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1915 if (nz>=fNdetectors || nz<0) {
1916 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1920 // ad hoc correction for 3rd ladder of SDD inner layer,
1921 // which is reversed (rotated by pi around local y)
1922 // this correction is OK only from AliITSv11Hybrid onwards
1923 if (GetR()>12. && GetR()<20.) { // SDD inner
1924 if(np==2) { // 3rd ladder
1925 Double_t posMod252[3];
1926 AliITSgeomTGeo::GetTranslation(252,posMod252);
1927 // check the Z coordinate of Mod 252: if negative
1928 // (old SDD geometry in AliITSv11Hybrid)
1929 // the swap of numeration whould be applied
1930 if(posMod252[2]<0.){
1931 nz = (fNdetectors-1) - nz;
1935 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1938 return np*fNdetectors + nz;
1940 //------------------------------------------------------------------------
1941 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1943 //--------------------------------------------------------------------
1944 // This function returns clusters within the "window"
1945 //--------------------------------------------------------------------
1947 if (fCurrentSlice<0) {
1948 Double_t rpi2 = 2.*fR*TMath::Pi();
1949 for (Int_t i=fI; i<fImax; i++) {
1952 if (fYmax<y) y -= rpi2;
1953 if (fYmin>y) y += rpi2;
1954 if (y<fYmin) continue;
1955 if (y>fYmax) continue;
1957 // skip clusters that are in "extended" road but they
1958 // 3sigma error does not touch the original road
1959 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
1960 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
1962 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
1965 return fClusters[i];
1968 for (Int_t i=fI; i<fImax; i++) {
1969 if (fYcs[i]<fYmin) continue;
1970 if (fYcs[i]>fYmax) continue;
1971 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
1972 ci=fClusterIndexCs[i];
1974 return fClustersCs[i];
1979 //------------------------------------------------------------------------
1980 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1982 //--------------------------------------------------------------------
1983 // This function returns the layer thickness at this point (units X0)
1984 //--------------------------------------------------------------------
1986 x0=AliITSRecoParam::GetX0Air();
1987 if (43<fR&&fR<45) { //SSD2
1990 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1991 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1992 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1993 for (Int_t i=0; i<12; i++) {
1994 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1995 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1999 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2000 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2004 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2005 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2008 if (37<fR&&fR<41) { //SSD1
2011 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2012 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2013 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2014 for (Int_t i=0; i<11; i++) {
2015 if (TMath::Abs(z-3.9*i)<0.15) {
2016 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2020 if (TMath::Abs(z+3.9*i)<0.15) {
2021 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2025 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2026 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2029 if (13<fR&&fR<26) { //SDD
2032 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2034 if (TMath::Abs(y-1.80)<0.55) {
2036 for (Int_t j=0; j<20; j++) {
2037 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2038 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2041 if (TMath::Abs(y+1.80)<0.55) {
2043 for (Int_t j=0; j<20; j++) {
2044 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2045 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2049 for (Int_t i=0; i<4; i++) {
2050 if (TMath::Abs(z-7.3*i)<0.60) {
2052 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2055 if (TMath::Abs(z+7.3*i)<0.60) {
2057 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2062 if (6<fR&&fR<8) { //SPD2
2063 Double_t dd=0.0063; x0=21.5;
2065 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2066 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2068 if (3<fR&&fR<5) { //SPD1
2069 Double_t dd=0.0063; x0=21.5;
2071 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2072 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2077 //------------------------------------------------------------------------
2078 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2080 fRmisal(det.fRmisal),
2082 fSinPhi(det.fSinPhi),
2083 fCosPhi(det.fCosPhi),
2089 fNChips(det.fNChips),
2090 fChipIsBad(det.fChipIsBad)
2094 //------------------------------------------------------------------------
2095 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2096 const AliITSDetTypeRec *detTypeRec)
2098 //--------------------------------------------------------------------
2099 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2100 //--------------------------------------------------------------------
2102 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2103 // while in the tracker they start from 0 for each layer
2104 for(Int_t il=0; il<ilayer; il++)
2105 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2108 if (ilayer==0 || ilayer==1) { // ---------- SPD
2110 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2112 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2115 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2119 // Get calibration from AliITSDetTypeRec
2120 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2121 calib->SetModuleIndex(idet);
2122 AliITSCalibration *calibSPDdead = 0;
2123 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2124 if (calib->IsBad() ||
2125 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2128 // printf("lay %d bad %d\n",ilayer,idet);
2131 // Get segmentation from AliITSDetTypeRec
2132 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2134 // Read info about bad chips
2135 fNChips = segm->GetMaximumChipIndex()+1;
2136 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2137 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2138 fChipIsBad = new Bool_t[fNChips];
2139 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2140 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2141 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2142 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2147 //------------------------------------------------------------------------
2148 Double_t AliITStrackerMI::GetEffectiveThickness()
2150 //--------------------------------------------------------------------
2151 // Returns the thickness between the current layer and the vertex (units X0)
2152 //--------------------------------------------------------------------
2155 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2156 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2157 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2161 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2162 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2166 Double_t xn=fgLayers[fI].GetR();
2167 for (Int_t i=0; i<fI; i++) {
2168 Double_t xi=fgLayers[i].GetR();
2169 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2175 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2176 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2179 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2180 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2184 //------------------------------------------------------------------------
2185 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2186 //-------------------------------------------------------------------
2187 // This function returns number of clusters within the "window"
2188 //--------------------------------------------------------------------
2190 for (Int_t i=fI; i<fN; i++) {
2191 const AliITSRecPoint *c=fClusters[i];
2192 if (c->GetZ() > fZmax) break;
2193 if (c->IsUsed()) continue;
2194 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2195 Double_t y=fR*det.GetPhi() + c->GetY();
2197 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2198 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2200 if (y<fYmin) continue;
2201 if (y>fYmax) continue;
2206 //------------------------------------------------------------------------
2207 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2208 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2210 //--------------------------------------------------------------------
2211 // This function refits the track "track" at the position "x" using
2212 // the clusters from "clusters"
2213 // If "extra"==kTRUE,
2214 // the clusters from overlapped modules get attached to "track"
2215 // If "planeff"==kTRUE,
2216 // special approach for plane efficiency evaluation is applyed
2217 //--------------------------------------------------------------------
2219 Int_t index[AliITSgeomTGeo::kNLayers];
2221 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2222 Int_t nc=clusters->GetNumberOfClusters();
2223 for (k=0; k<nc; k++) {
2224 Int_t idx=clusters->GetClusterIndex(k);
2225 Int_t ilayer=(idx&0xf0000000)>>28;
2229 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2231 //------------------------------------------------------------------------
2232 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2233 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2235 //--------------------------------------------------------------------
2236 // This function refits the track "track" at the position "x" using
2237 // the clusters from array
2238 // If "extra"==kTRUE,
2239 // the clusters from overlapped modules get attached to "track"
2240 // If "planeff"==kTRUE,
2241 // special approach for plane efficiency evaluation is applyed
2242 //--------------------------------------------------------------------
2243 Int_t index[AliITSgeomTGeo::kNLayers];
2245 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2247 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2248 index[k]=clusters[k];
2251 // special for cosmics: check which the innermost layer crossed
2253 Int_t innermostlayer=5;
2254 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2255 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2256 if(drphi < fgLayers[innermostlayer].GetR()) break;
2258 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2260 Int_t modstatus=1; // found
2262 Int_t from, to, step;
2263 if (xx > track->GetX()) {
2264 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2267 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2270 TString dir = (step>0 ? "outward" : "inward");
2272 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2273 AliITSlayer &layer=fgLayers[ilayer];
2274 Double_t r=layer.GetR();
2276 if (step<0 && xx>r) break;
2278 // material between SSD and SDD, SDD and SPD
2279 Double_t hI=ilayer-0.5*step;
2280 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2281 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2282 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2283 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2286 Double_t oldGlobXYZ[3];
2287 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2289 // continue if we are already beyond this layer
2290 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2291 if(step>0 && oldGlobR > r) continue; // going outward
2292 if(step<0 && oldGlobR < r) continue; // going inward
2295 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2297 Int_t idet=layer.FindDetectorIndex(phi,z);
2299 // check if we allow a prolongation without point for large-eta tracks
2300 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2302 modstatus = 4; // out in z
2303 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2304 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2307 // apply correction for material of the current layer
2308 // add time if going outward
2309 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2313 if (idet<0) return kFALSE;
2316 const AliITSdetector &det=layer.GetDetector(idet);
2317 // only for ITS-SA tracks refit
2318 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2320 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2322 track->SetDetectorIndex(idet);
2323 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2325 Double_t dz,zmin,zmax,dy,ymin,ymax;
2327 const AliITSRecPoint *clAcc=0;
2328 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2330 Int_t idx=index[ilayer];
2331 if (idx>=0) { // cluster in this layer
2332 modstatus = 6; // no refit
2333 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2335 if (idet != cl->GetDetectorIndex()) {
2336 idet=cl->GetDetectorIndex();
2337 const AliITSdetector &detc=layer.GetDetector(idet);
2338 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2339 track->SetDetectorIndex(idet);
2340 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2342 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2343 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2347 modstatus = 1; // found
2352 } else { // no cluster in this layer
2354 modstatus = 3; // skipped
2355 // Plane Eff determination:
2356 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2357 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2358 UseTrackForPlaneEff(track,ilayer);
2361 modstatus = 5; // no cls in road
2363 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2364 dz = 0.5*(zmax-zmin);
2365 dy = 0.5*(ymax-ymin);
2366 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2367 if (dead==1) modstatus = 7; // holes in z in SPD
2368 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2373 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2374 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2376 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2379 if (extra) { // search for extra clusters in overlapped modules
2380 AliITStrackV2 tmp(*track);
2381 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2382 layer.SelectClusters(zmin,zmax,ymin,ymax);
2384 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2386 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2387 Double_t tolerance=0.1;
2388 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2389 // only clusters in another module! (overlaps)
2390 idetExtra = clExtra->GetDetectorIndex();
2391 if (idet == idetExtra) continue;
2393 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2395 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2396 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2397 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2398 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2400 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2401 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2404 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2405 track->SetExtraModule(ilayer,idetExtra);
2407 } // end search for extra clusters in overlapped modules
2409 // Correct for material of the current layer
2411 // add time if going outward
2412 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2413 track->SetCheckInvariant(kTRUE);
2414 } // end loop on layers
2416 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2420 //------------------------------------------------------------------------
2421 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2424 // calculate normalized chi2
2425 // return NormalizedChi2(track,0);
2428 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2429 // track->fdEdxMismatch=0;
2430 Float_t dedxmismatch =0;
2431 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2433 for (Int_t i = 0;i<6;i++){
2434 if (track->GetClIndex(i)>=0){
2435 Float_t cerry, cerrz;
2436 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2438 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2441 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2442 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2443 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2445 cchi2+=(0.5-ratio)*10.;
2446 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2447 dedxmismatch+=(0.5-ratio)*10.;
2451 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2452 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2453 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2454 if (i<2) chi2+=2*cl->GetDeltaProbability();
2460 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2461 track->SetdEdxMismatch(dedxmismatch);
2465 for (Int_t i = 0;i<4;i++){
2466 if (track->GetClIndex(i)>=0){
2467 Float_t cerry, cerrz;
2468 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2469 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2472 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2473 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2477 for (Int_t i = 4;i<6;i++){
2478 if (track->GetClIndex(i)>=0){
2479 Float_t cerry, cerrz;
2480 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2481 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2484 Float_t cerryb, cerrzb;
2485 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2486 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2489 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2490 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2495 if (track->GetESDtrack()->GetTPCsignal()>85){
2496 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2498 chi2+=(0.5-ratio)*5.;
2501 chi2+=(ratio-2.0)*3;
2505 Double_t match = TMath::Sqrt(track->GetChi22());
2506 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2507 if (!track->GetConstrain()) {
2508 if (track->GetNumberOfClusters()>2) {
2509 match/=track->GetNumberOfClusters()-2.;
2514 if (match<0) match=0;
2516 // penalty factor for missing points (NDeadZone>0), but no penalty
2517 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2518 // or there is a dead from OCDB)
2519 Float_t deadzonefactor = 0.;
2520 if (track->GetNDeadZone()>0.) {
2521 Int_t sumDeadZoneProbability=0;
2522 for(Int_t ilay=0;ilay<6;ilay++) {
2523 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2525 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2526 if(nDeadZoneWithProbNot1>0) {
2527 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2528 AliDebug(2,Form("nDeadZone %f sumDZProbability %f nDZWithProbNot1 %f deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2529 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2531 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2532 deadzonefactor = 3.*(1.1-deadZoneProbability);
2536 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2537 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2538 1./(1.+track->GetNSkipped()));
2539 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2540 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2543 //------------------------------------------------------------------------
2544 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2547 // return matching chi2 between two tracks
2548 Double_t largeChi2=1000.;
2550 AliITStrackMI track3(*track2);
2551 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2553 vec(0,0)=track1->GetY() - track3.GetY();
2554 vec(1,0)=track1->GetZ() - track3.GetZ();
2555 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2556 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2557 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2560 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2561 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2562 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2563 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2564 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2566 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2567 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2568 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2569 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2571 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2572 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2573 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2575 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2576 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2578 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2581 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2582 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2585 //------------------------------------------------------------------------
2586 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2589 // return probability that given point (characterized by z position and error)
2590 // is in SPD dead zone
2591 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2593 Double_t probability = 0.;
2594 Double_t nearestz = 0.,distz=0.;
2595 Int_t nearestzone = -1;
2596 Double_t mindistz = 1000.;
2598 // find closest dead zone
2599 for (Int_t i=0; i<3; i++) {
2600 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2601 if (distz<mindistz) {
2603 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2608 // too far from dead zone
2609 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2612 Double_t zmin, zmax;
2613 if (nearestzone==0) { // dead zone at z = -7
2614 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2615 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2616 } else if (nearestzone==1) { // dead zone at z = 0
2617 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2618 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2619 } else if (nearestzone==2) { // dead zone at z = +7
2620 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2621 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2626 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2628 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2629 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2630 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2633 //------------------------------------------------------------------------
2634 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2637 // calculate normalized chi2
2639 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2641 for (Int_t i = 0;i<6;i++){
2642 if (TMath::Abs(track->GetDy(i))>0){
2643 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2644 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2647 else{chi2[i]=10000;}
2650 TMath::Sort(6,chi2,index,kFALSE);
2651 Float_t max = float(ncl)*fac-1.;
2652 Float_t sumchi=0, sumweight=0;
2653 for (Int_t i=0;i<max+1;i++){
2654 Float_t weight = (i<max)?1.:(max+1.-i);
2655 sumchi+=weight*chi2[index[i]];
2658 Double_t normchi2 = sumchi/sumweight;
2661 //------------------------------------------------------------------------
2662 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2665 // calculate normalized chi2
2666 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2669 for (Int_t i=0;i<6;i++){
2670 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2671 Double_t sy1 = forwardtrack->GetSigmaY(i);
2672 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2673 Double_t sy2 = backtrack->GetSigmaY(i);
2674 Double_t sz2 = backtrack->GetSigmaZ(i);
2675 if (i<2){ sy2=1000.;sz2=1000;}
2677 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2678 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2680 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2681 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2683 res+= nz0*nz0+ny0*ny0;
2686 if (npoints>1) return
2687 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2688 //2*forwardtrack->fNUsed+
2689 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2690 1./(1.+forwardtrack->GetNSkipped()));
2693 //------------------------------------------------------------------------
2694 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2695 //--------------------------------------------------------------------
2696 // Return pointer to a given cluster
2697 //--------------------------------------------------------------------
2698 Int_t l=(index & 0xf0000000) >> 28;
2699 Int_t c=(index & 0x0fffffff) >> 00;
2700 return fgLayers[l].GetWeight(c);
2702 //------------------------------------------------------------------------
2703 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2705 //---------------------------------------------
2706 // register track to the list
2708 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2711 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2712 Int_t index = track->GetClusterIndex(icluster);
2713 Int_t l=(index & 0xf0000000) >> 28;
2714 Int_t c=(index & 0x0fffffff) >> 00;
2715 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2716 for (Int_t itrack=0;itrack<4;itrack++){
2717 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2718 fgLayers[l].SetClusterTracks(itrack,c,id);
2724 //------------------------------------------------------------------------
2725 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2727 //---------------------------------------------
2728 // unregister track from the list
2729 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2730 Int_t index = track->GetClusterIndex(icluster);
2731 Int_t l=(index & 0xf0000000) >> 28;
2732 Int_t c=(index & 0x0fffffff) >> 00;
2733 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2734 for (Int_t itrack=0;itrack<4;itrack++){
2735 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2736 fgLayers[l].SetClusterTracks(itrack,c,-1);
2741 //------------------------------------------------------------------------
2742 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2744 //-------------------------------------------------------------
2745 //get number of shared clusters
2746 //-------------------------------------------------------------
2748 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2749 // mean number of clusters
2750 Float_t *ny = GetNy(id), *nz = GetNz(id);
2753 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2754 Int_t index = track->GetClusterIndex(icluster);
2755 Int_t l=(index & 0xf0000000) >> 28;
2756 Int_t c=(index & 0x0fffffff) >> 00;
2757 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2759 printf("problem\n");
2761 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2765 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2766 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2767 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2769 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2772 deltan = (cl->GetNz()-nz[l]);
2774 if (deltan>2.0) continue; // extended - highly probable shared cluster
2775 weight = 2./TMath::Max(3.+deltan,2.);
2777 for (Int_t itrack=0;itrack<4;itrack++){
2778 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2780 clist[l] = (AliITSRecPoint*)GetCluster(index);
2786 track->SetNUsed(shared);
2789 //------------------------------------------------------------------------
2790 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2793 // find first shared track
2795 // mean number of clusters
2796 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2798 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2799 Int_t sharedtrack=100000;
2800 Int_t tracks[24],trackindex=0;
2801 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2803 for (Int_t icluster=0;icluster<6;icluster++){
2804 if (clusterlist[icluster]<0) continue;
2805 Int_t index = clusterlist[icluster];
2806 Int_t l=(index & 0xf0000000) >> 28;
2807 Int_t c=(index & 0x0fffffff) >> 00;
2809 printf("problem\n");
2811 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2812 //if (l>3) continue;
2813 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2816 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2817 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2818 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2820 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2823 deltan = (cl->GetNz()-nz[l]);
2825 if (deltan>2.0) continue; // extended - highly probable shared cluster
2827 for (Int_t itrack=3;itrack>=0;itrack--){
2828 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2829 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2830 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2835 if (trackindex==0) return -1;
2837 sharedtrack = tracks[0];
2839 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2842 Int_t tracks2[24], cluster[24];
2843 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2846 for (Int_t i=0;i<trackindex;i++){
2847 if (tracks[i]<0) continue;
2848 tracks2[index] = tracks[i];
2850 for (Int_t j=i+1;j<trackindex;j++){
2851 if (tracks[j]<0) continue;
2852 if (tracks[j]==tracks[i]){
2860 for (Int_t i=0;i<index;i++){
2861 if (cluster[index]>max) {
2862 sharedtrack=tracks2[index];
2869 if (sharedtrack>=100000) return -1;
2871 // make list of overlaps
2873 for (Int_t icluster=0;icluster<6;icluster++){
2874 if (clusterlist[icluster]<0) continue;
2875 Int_t index = clusterlist[icluster];
2876 Int_t l=(index & 0xf0000000) >> 28;
2877 Int_t c=(index & 0x0fffffff) >> 00;
2878 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2879 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2881 if (cl->GetNy()>2) continue;
2882 if (cl->GetNz()>2) continue;
2885 if (cl->GetNy()>3) continue;
2886 if (cl->GetNz()>3) continue;
2889 for (Int_t itrack=3;itrack>=0;itrack--){
2890 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2891 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2899 //------------------------------------------------------------------------
2900 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2902 // try to find track hypothesys without conflicts
2903 // with minimal chi2;
2904 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2905 Int_t entries1 = arr1->GetEntriesFast();
2906 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2907 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2908 Int_t entries2 = arr2->GetEntriesFast();
2909 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2911 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2912 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2913 if (track10->Pt()>0.5+track20->Pt()) return track10;
2915 for (Int_t itrack=0;itrack<entries1;itrack++){
2916 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2917 UnRegisterClusterTracks(track,trackID1);
2920 for (Int_t itrack=0;itrack<entries2;itrack++){
2921 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2922 UnRegisterClusterTracks(track,trackID2);
2926 Float_t maxconflicts=6;
2927 Double_t maxchi2 =1000.;
2929 // get weight of hypothesys - using best hypothesy
2932 Int_t list1[6],list2[6];
2933 AliITSRecPoint *clist1[6], *clist2[6] ;
2934 RegisterClusterTracks(track10,trackID1);
2935 RegisterClusterTracks(track20,trackID2);
2936 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2937 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2938 UnRegisterClusterTracks(track10,trackID1);
2939 UnRegisterClusterTracks(track20,trackID2);
2942 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2943 Float_t nerry[6],nerrz[6];
2944 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2945 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2946 for (Int_t i=0;i<6;i++){
2947 if ( (erry1[i]>0) && (erry2[i]>0)) {
2948 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2949 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2951 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2952 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2954 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2955 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2956 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2959 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2960 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2961 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2969 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2970 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2971 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2972 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2974 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2975 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2976 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2978 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2979 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2980 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2983 Double_t sumw = w1+w2;
2987 w1 = (d2+0.5)/(d1+d2+1.);
2988 w2 = (d1+0.5)/(d1+d2+1.);
2990 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2991 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2993 // get pair of "best" hypothesys
2995 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2996 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2998 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2999 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3000 //if (track1->fFakeRatio>0) continue;
3001 RegisterClusterTracks(track1,trackID1);
3002 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3003 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3005 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3006 //if (track2->fFakeRatio>0) continue;
3008 RegisterClusterTracks(track2,trackID2);
3009 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3010 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3011 UnRegisterClusterTracks(track2,trackID2);
3013 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3014 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3015 if (nskipped>0.5) continue;
3017 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3018 if (conflict1+1<cconflict1) continue;
3019 if (conflict2+1<cconflict2) continue;
3023 for (Int_t i=0;i<6;i++){
3025 Float_t c1 =0.; // conflict coeficients
3027 if (clist1[i]&&clist2[i]){
3030 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3033 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3035 c1 = 2./TMath::Max(3.+deltan,2.);
3036 c2 = 2./TMath::Max(3.+deltan,2.);
3042 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3045 deltan = (clist1[i]->GetNz()-nz1[i]);
3047 c1 = 2./TMath::Max(3.+deltan,2.);
3054 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3057 deltan = (clist2[i]->GetNz()-nz2[i]);
3059 c2 = 2./TMath::Max(3.+deltan,2.);
3065 if (TMath::Abs(track1->GetDy(i))>0.) {
3066 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3067 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3068 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3069 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3071 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3074 if (TMath::Abs(track2->GetDy(i))>0.) {
3075 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3076 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3077 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3078 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3081 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3083 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3084 if (chi21>0) sum+=w1;
3085 if (chi22>0) sum+=w2;
3088 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3089 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3090 Double_t normchi2 = 2*conflict+sumchi2/sum;
3091 if ( normchi2 <maxchi2 ){
3094 maxconflicts = conflict;
3098 UnRegisterClusterTracks(track1,trackID1);
3101 // if (maxconflicts<4 && maxchi2<th0){
3102 if (maxchi2<th0*2.){
3103 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3104 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3105 track1->SetChi2MIP(5,maxconflicts);
3106 track1->SetChi2MIP(6,maxchi2);
3107 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3108 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3109 track1->SetChi2MIP(8,index1);
3110 fBestTrackIndex[trackID1] =index1;
3111 UpdateESDtrack(track1, AliESDtrack::kITSin);
3113 else if (track10->GetChi2MIP(0)<th1){
3114 track10->SetChi2MIP(5,maxconflicts);
3115 track10->SetChi2MIP(6,maxchi2);
3116 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3117 UpdateESDtrack(track10,AliESDtrack::kITSin);
3120 for (Int_t itrack=0;itrack<entries1;itrack++){
3121 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3122 UnRegisterClusterTracks(track,trackID1);
3125 for (Int_t itrack=0;itrack<entries2;itrack++){
3126 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3127 UnRegisterClusterTracks(track,trackID2);
3130 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3131 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3132 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3133 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3134 RegisterClusterTracks(track10,trackID1);
3136 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3137 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3138 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3139 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3140 RegisterClusterTracks(track20,trackID2);
3145 //------------------------------------------------------------------------
3146 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3147 //--------------------------------------------------------------------
3148 // This function marks clusters assigned to the track
3149 //--------------------------------------------------------------------
3150 AliTracker::UseClusters(t,from);
3152 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3153 //if (c->GetQ()>2) c->Use();
3154 if (c->GetSigmaZ2()>0.1) c->Use();
3155 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3156 //if (c->GetQ()>2) c->Use();
3157 if (c->GetSigmaZ2()>0.1) c->Use();
3160 //------------------------------------------------------------------------
3161 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3163 //------------------------------------------------------------------
3164 // add track to the list of hypothesys
3165 //------------------------------------------------------------------
3167 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3168 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3170 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3172 array = new TObjArray(10);
3173 fTrackHypothesys.AddAt(array,esdindex);
3175 array->AddLast(track);
3177 //------------------------------------------------------------------------
3178 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3180 //-------------------------------------------------------------------
3181 // compress array of track hypothesys
3182 // keep only maxsize best hypothesys
3183 //-------------------------------------------------------------------
3184 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3185 if (! (fTrackHypothesys.At(esdindex)) ) return;
3186 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3187 Int_t entries = array->GetEntriesFast();
3189 //- find preliminary besttrack as a reference
3190 Float_t minchi2=10000;
3192 AliITStrackMI * besttrack=0;
3193 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3194 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3195 if (!track) continue;
3196 Float_t chi2 = NormalizedChi2(track,0);
3198 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3199 track->SetLabel(tpcLabel);
3201 track->SetFakeRatio(1.);
3202 CookLabel(track,0.); //For comparison only
3204 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3205 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3206 if (track->GetNumberOfClusters()<maxn) continue;
3207 maxn = track->GetNumberOfClusters();
3214 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3215 delete array->RemoveAt(itrack);
3219 if (!besttrack) return;
3222 //take errors of best track as a reference
3223 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3224 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3225 for (Int_t j=0;j<6;j++) {
3226 if (besttrack->GetClIndex(j)>=0){
3227 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3228 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3229 ny[j] = besttrack->GetNy(j);
3230 nz[j] = besttrack->GetNz(j);
3234 // calculate normalized chi2
3236 Float_t * chi2 = new Float_t[entries];
3237 Int_t * index = new Int_t[entries];
3238 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3239 for (Int_t itrack=0;itrack<entries;itrack++){
3240 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3242 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3243 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3244 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3245 chi2[itrack] = track->GetChi2MIP(0);
3247 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3248 delete array->RemoveAt(itrack);
3254 TMath::Sort(entries,chi2,index,kFALSE);
3255 besttrack = (AliITStrackMI*)array->At(index[0]);
3256 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3257 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3258 for (Int_t j=0;j<6;j++){
3259 if (besttrack->GetClIndex(j)>=0){
3260 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3261 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3262 ny[j] = besttrack->GetNy(j);
3263 nz[j] = besttrack->GetNz(j);
3268 // calculate one more time with updated normalized errors
3269 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3270 for (Int_t itrack=0;itrack<entries;itrack++){
3271 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3273 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3274 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3275 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3276 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3279 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3280 delete array->RemoveAt(itrack);
3285 entries = array->GetEntriesFast();
3289 TObjArray * newarray = new TObjArray();
3290 TMath::Sort(entries,chi2,index,kFALSE);
3291 besttrack = (AliITStrackMI*)array->At(index[0]);
3293 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3295 for (Int_t j=0;j<6;j++){
3296 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3297 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3298 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3299 ny[j] = besttrack->GetNy(j);
3300 nz[j] = besttrack->GetNz(j);
3303 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3304 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3305 Float_t minn = besttrack->GetNumberOfClusters()-3;
3307 for (Int_t i=0;i<entries;i++){
3308 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3309 if (!track) continue;
3310 if (accepted>maxcut) break;
3311 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3312 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3313 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3314 delete array->RemoveAt(index[i]);
3318 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3319 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3320 if (!shortbest) accepted++;
3322 newarray->AddLast(array->RemoveAt(index[i]));
3323 for (Int_t j=0;j<6;j++){
3325 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3326 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3327 ny[j] = track->GetNy(j);
3328 nz[j] = track->GetNz(j);
3333 delete array->RemoveAt(index[i]);
3337 delete fTrackHypothesys.RemoveAt(esdindex);
3338 fTrackHypothesys.AddAt(newarray,esdindex);
3342 delete fTrackHypothesys.RemoveAt(esdindex);
3348 //------------------------------------------------------------------------
3349 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3351 //-------------------------------------------------------------
3352 // try to find best hypothesy
3353 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3354 //-------------------------------------------------------------
3355 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3356 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3357 if (!array) return 0;
3358 Int_t entries = array->GetEntriesFast();
3359 if (!entries) return 0;
3360 Float_t minchi2 = 100000;
3361 AliITStrackMI * besttrack=0;
3363 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3364 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3365 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3366 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3368 for (Int_t i=0;i<entries;i++){
3369 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3370 if (!track) continue;
3371 Float_t sigmarfi,sigmaz;
3372 GetDCASigma(track,sigmarfi,sigmaz);
3373 track->SetDnorm(0,sigmarfi);
3374 track->SetDnorm(1,sigmaz);
3376 track->SetChi2MIP(1,1000000);
3377 track->SetChi2MIP(2,1000000);
3378 track->SetChi2MIP(3,1000000);
3381 backtrack = new(backtrack) AliITStrackMI(*track);
3382 if (track->GetConstrain()) {
3383 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3384 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3385 backtrack->ResetCovariance(10.);
3387 backtrack->ResetCovariance(10.);
3389 backtrack->ResetClusters();
3391 Double_t x = original->GetX();
3392 if (!RefitAt(x,backtrack,track)) continue;
3394 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3395 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3396 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3397 track->SetChi22(GetMatchingChi2(backtrack,original));
3399 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3400 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3401 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3404 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3406 //forward track - without constraint
3407 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3408 forwardtrack->ResetClusters();
3410 RefitAt(x,forwardtrack,track);
3411 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3412 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3413 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3415 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3416 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3417 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3418 forwardtrack->SetD(0,track->GetD(0));
3419 forwardtrack->SetD(1,track->GetD(1));
3422 AliITSRecPoint* clist[6];
3423 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3424 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3427 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3428 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3429 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3430 track->SetChi2MIP(3,1000);
3433 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3435 for (Int_t ichi=0;ichi<5;ichi++){
3436 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3438 if (chi2 < minchi2){
3439 //besttrack = new AliITStrackMI(*forwardtrack);
3441 besttrack->SetLabel(track->GetLabel());
3442 besttrack->SetFakeRatio(track->GetFakeRatio());
3444 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3445 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3446 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3450 delete forwardtrack;
3452 for (Int_t i=0;i<entries;i++){
3453 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3455 if (!track) continue;
3457 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3458 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3459 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3460 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3461 delete array->RemoveAt(i);
3471 SortTrackHypothesys(esdindex,checkmax,1);
3473 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3474 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3475 besttrack = (AliITStrackMI*)array->At(0);
3476 if (!besttrack) return 0;
3477 besttrack->SetChi2MIP(8,0);
3478 fBestTrackIndex[esdindex]=0;
3479 entries = array->GetEntriesFast();
3480 AliITStrackMI *longtrack =0;
3482 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3483 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3484 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3485 if (!track->GetConstrain()) continue;
3486 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3487 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3488 if (track->GetChi2MIP(0)>4.) continue;
3489 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3492 //if (longtrack) besttrack=longtrack;
3495 AliITSRecPoint * clist[6];
3496 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3497 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3498 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3499 RegisterClusterTracks(besttrack,esdindex);
3506 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3507 if (sharedtrack>=0){
3509 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3511 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3517 if (shared>2.5) return 0;
3518 if (shared>1.0) return besttrack;
3520 // Don't sign clusters if not gold track
3522 if (!besttrack->IsGoldPrimary()) return besttrack;
3523 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3525 if (fConstraint[fPass]){
3529 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3530 for (Int_t i=0;i<6;i++){
3531 Int_t index = besttrack->GetClIndex(i);
3532 if (index<0) continue;
3533 Int_t ilayer = (index & 0xf0000000) >> 28;
3534 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3535 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3537 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3538 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3539 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3540 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3541 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3542 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3544 Bool_t cansign = kTRUE;
3545 for (Int_t itrack=0;itrack<entries; itrack++){
3546 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3547 if (!track) continue;
3548 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3549 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3555 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3556 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3557 if (!c->IsUsed()) c->Use();
3563 //------------------------------------------------------------------------
3564 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3567 // get "best" hypothesys
3570 Int_t nentries = itsTracks.GetEntriesFast();
3571 for (Int_t i=0;i<nentries;i++){
3572 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3573 if (!track) continue;
3574 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3575 if (!array) continue;
3576 if (array->GetEntriesFast()<=0) continue;
3578 AliITStrackMI* longtrack=0;
3580 Float_t maxchi2=1000;
3581 for (Int_t j=0;j<array->GetEntriesFast();j++){
3582 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3583 if (!trackHyp) continue;
3584 if (trackHyp->GetGoldV0()) {
3585 longtrack = trackHyp; //gold V0 track taken
3588 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3589 Float_t chi2 = trackHyp->GetChi2MIP(0);
3591 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3593 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3595 if (chi2 > maxchi2) continue;
3596 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3603 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3604 if (!longtrack) {longtrack = besttrack;}
3605 else besttrack= longtrack;
3609 AliITSRecPoint * clist[6];
3610 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3612 track->SetNUsed(shared);
3613 track->SetNSkipped(besttrack->GetNSkipped());
3614 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3616 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3620 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3621 //if (sharedtrack==-1) sharedtrack=0;
3622 if (sharedtrack>=0) {
3623 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3626 if (besttrack&&fAfterV0) {
3627 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3629 if (besttrack&&fConstraint[fPass])
3630 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3631 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3632 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3633 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3639 //------------------------------------------------------------------------
3640 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3641 //--------------------------------------------------------------------
3642 //This function "cooks" a track label. If label<0, this track is fake.
3643 //--------------------------------------------------------------------
3646 if (track->GetESDtrack()){
3647 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3648 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3649 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3651 track->SetChi2MIP(9,0);
3653 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3654 Int_t cindex = track->GetClusterIndex(i);
3655 Int_t l=(cindex & 0xf0000000) >> 28;
3656 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3658 for (Int_t ind=0;ind<3;ind++){
3659 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3660 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3662 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3665 Int_t nclusters = track->GetNumberOfClusters();
3666 if (nclusters > 0) //PH Some tracks don't have any cluster
3667 track->SetFakeRatio(double(nwrong)/double(nclusters));
3668 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3669 track->SetLabel(-tpcLabel);
3671 track->SetLabel(tpcLabel);
3673 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3676 //------------------------------------------------------------------------
3677 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
3679 // Fill the dE/dx in this track
3681 track->SetChi2MIP(9,0);
3682 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3683 Int_t cindex = track->GetClusterIndex(i);
3684 Int_t l=(cindex & 0xf0000000) >> 28;
3685 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3686 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3688 for (Int_t ind=0;ind<3;ind++){
3689 if (cl->GetLabel(ind)==lab) isWrong=0;
3691 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3695 track->CookdEdx(low,up);
3697 //------------------------------------------------------------------------
3698 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3700 // Create some arrays
3702 if (fCoefficients) delete []fCoefficients;
3703 fCoefficients = new Float_t[ntracks*48];
3704 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3706 //------------------------------------------------------------------------
3707 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3710 // Compute predicted chi2
3712 Float_t erry,errz,covyz;
3713 Float_t theta = track->GetTgl();
3714 Float_t phi = track->GetSnp();
3715 phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3716 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3717 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()));
3718 // Take into account the mis-alignment (bring track to cluster plane)
3719 Double_t xTrOrig=track->GetX();
3720 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3721 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()));
3722 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3723 // Bring the track back to detector plane in ideal geometry
3724 // [mis-alignment will be accounted for in UpdateMI()]
3725 if (!track->Propagate(xTrOrig)) return 1000.;
3727 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3728 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3730 chi2+=0.5*TMath::Min(delta/2,2.);
3731 chi2+=2.*cluster->GetDeltaProbability();
3734 track->SetNy(layer,ny);
3735 track->SetNz(layer,nz);
3736 track->SetSigmaY(layer,erry);
3737 track->SetSigmaZ(layer, errz);
3738 track->SetSigmaYZ(layer,covyz);
3739 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3740 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3744 //------------------------------------------------------------------------
3745 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3750 Int_t layer = (index & 0xf0000000) >> 28;
3751 track->SetClIndex(layer, index);
3752 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3753 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3754 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3755 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3759 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
3762 // Take into account the mis-alignment (bring track to cluster plane)
3763 Double_t xTrOrig=track->GetX();
3764 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3765 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3766 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3767 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3769 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3772 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3773 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3774 c.SetSigmaYZ(track->GetSigmaYZ(layer));
3777 Int_t updated = track->UpdateMI(&c,chi2,index);
3778 // Bring the track back to detector plane in ideal geometry
3779 if (!track->Propagate(xTrOrig)) return 0;
3781 if(!updated) AliDebug(2,"update failed");
3785 //------------------------------------------------------------------------
3786 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3789 //DCA sigmas parameterization
3790 //to be paramterized using external parameters in future
3793 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3794 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3796 //------------------------------------------------------------------------
3797 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3800 // Clusters from delta electrons?
3802 Int_t entries = clusterArray->GetEntriesFast();
3803 if (entries<4) return;
3804 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3805 Int_t layer = cluster->GetLayer();
3806 if (layer>1) return;
3808 Int_t ncandidates=0;
3809 Float_t r = (layer>0)? 7:4;
3811 for (Int_t i=0;i<entries;i++){
3812 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3813 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3814 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3815 index[ncandidates] = i; //candidate to belong to delta electron track
3817 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3818 cl0->SetDeltaProbability(1);
3824 for (Int_t i=0;i<ncandidates;i++){
3825 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3826 if (cl0->GetDeltaProbability()>0.8) continue;
3829 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3830 sumy=sumz=sumy2=sumyz=sumw=0.0;
3831 for (Int_t j=0;j<ncandidates;j++){
3833 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3835 Float_t dz = cl0->GetZ()-cl1->GetZ();
3836 Float_t dy = cl0->GetY()-cl1->GetY();
3837 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3838 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3839 y[ncl] = cl1->GetY();
3840 z[ncl] = cl1->GetZ();
3841 sumy+= y[ncl]*weight;
3842 sumz+= z[ncl]*weight;
3843 sumy2+=y[ncl]*y[ncl]*weight;
3844 sumyz+=y[ncl]*z[ncl]*weight;
3849 if (ncl<4) continue;
3850 Float_t det = sumw*sumy2 - sumy*sumy;
3852 if (TMath::Abs(det)>0.01){
3853 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3854 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3855 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3858 Float_t z0 = sumyz/sumy;
3859 delta = TMath::Abs(cl0->GetZ()-z0);
3862 cl0->SetDeltaProbability(1-20.*delta);
3866 //------------------------------------------------------------------------
3867 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3872 track->UpdateESDtrack(flags);
3873 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3874 if (oldtrack) delete oldtrack;
3875 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3876 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3877 printf("Problem\n");
3880 //------------------------------------------------------------------------
3881 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3883 // Get nearest upper layer close to the point xr.
3884 // rough approximation
3886 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3887 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3889 for (Int_t i=0;i<6;i++){
3890 if (radius<kRadiuses[i]){
3897 //------------------------------------------------------------------------
3898 void AliITStrackerMI::BuildMaterialLUT(TString material) {
3899 //--------------------------------------------------------------------
3900 // Fill a look-up table with mean material
3901 //--------------------------------------------------------------------
3905 Double_t point1[3],point2[3];
3906 Double_t phi,cosphi,sinphi,z;
3907 // 0-5 layers, 6 pipe, 7-8 shields
3908 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
3909 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
3911 Int_t ifirst=0,ilast=0;
3912 if(material.Contains("Pipe")) {
3914 } else if(material.Contains("Shields")) {
3916 } else if(material.Contains("Layers")) {
3919 Error("BuildMaterialLUT","Wrong layer name\n");
3922 for(Int_t imat=ifirst; imat<=ilast; imat++) {
3923 Double_t param[5]={0.,0.,0.,0.,0.};
3924 for (Int_t i=0; i<n; i++) {
3925 phi = 2.*TMath::Pi()*gRandom->Rndm();
3926 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
3927 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
3928 point1[0] = rmin[imat]*cosphi;
3929 point1[1] = rmin[imat]*sinphi;
3931 point2[0] = rmax[imat]*cosphi;
3932 point2[1] = rmax[imat]*sinphi;
3934 AliTracker::MeanMaterialBudget(point1,point2,mparam);
3935 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
3937 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
3939 fxOverX0Layer[imat] = param[1];
3940 fxTimesRhoLayer[imat] = param[0]*param[4];
3941 } else if(imat==6) {
3942 fxOverX0Pipe = param[1];
3943 fxTimesRhoPipe = param[0]*param[4];
3944 } else if(imat==7) {
3945 fxOverX0Shield[0] = param[1];
3946 fxTimesRhoShield[0] = param[0]*param[4];
3947 } else if(imat==8) {
3948 fxOverX0Shield[1] = param[1];
3949 fxTimesRhoShield[1] = param[0]*param[4];
3953 printf("%s\n",material.Data());
3954 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
3955 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
3956 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
3957 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
3958 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
3959 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
3960 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
3961 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
3962 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
3966 //------------------------------------------------------------------------
3967 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
3968 TString direction) {
3969 //-------------------------------------------------------------------
3970 // Propagate beyond beam pipe and correct for material
3971 // (material budget in different ways according to fUseTGeo value)
3972 // Add time if going outward (PropagateTo or PropagateToTGeo)
3973 //-------------------------------------------------------------------
3975 // Define budget mode:
3976 // 0: material from AliITSRecoParam (hard coded)
3977 // 1: material from TGeo in one step (on the fly)
3978 // 2: material from lut
3979 // 3: material from TGeo in one step (same for all hypotheses)
3992 if(fTrackingPhase.Contains("Clusters2Tracks"))
3993 { mode=3; } else { mode=1; }
3996 if(fTrackingPhase.Contains("Clusters2Tracks"))
3997 { mode=3; } else { mode=2; }
4003 if(fTrackingPhase.Contains("Default")) mode=0;
4005 Int_t index=fCurrentEsdTrack;
4007 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4008 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4010 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4012 Double_t xOverX0,x0,lengthTimesMeanDensity;
4016 xOverX0 = AliITSRecoParam::GetdPipe();
4017 x0 = AliITSRecoParam::GetX0Be();
4018 lengthTimesMeanDensity = xOverX0*x0;
4019 lengthTimesMeanDensity *= dir;
4020 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4023 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4026 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4027 xOverX0 = fxOverX0Pipe;
4028 lengthTimesMeanDensity = fxTimesRhoPipe;
4029 lengthTimesMeanDensity *= dir;
4030 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4033 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4034 if(fxOverX0PipeTrks[index]<0) {
4035 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4036 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4037 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4038 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4039 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4042 xOverX0 = fxOverX0PipeTrks[index];
4043 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4044 lengthTimesMeanDensity *= dir;
4045 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4051 //------------------------------------------------------------------------
4052 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4054 TString direction) {
4055 //-------------------------------------------------------------------
4056 // Propagate beyond SPD or SDD shield and correct for material
4057 // (material budget in different ways according to fUseTGeo value)
4058 // Add time if going outward (PropagateTo or PropagateToTGeo)
4059 //-------------------------------------------------------------------
4061 // Define budget mode:
4062 // 0: material from AliITSRecoParam (hard coded)
4063 // 1: material from TGeo in steps of X cm (on the fly)
4064 // X = AliITSRecoParam::GetStepSizeTGeo()
4065 // 2: material from lut
4066 // 3: material from TGeo in one step (same for all hypotheses)
4079 if(fTrackingPhase.Contains("Clusters2Tracks"))
4080 { mode=3; } else { mode=1; }
4083 if(fTrackingPhase.Contains("Clusters2Tracks"))
4084 { mode=3; } else { mode=2; }
4090 if(fTrackingPhase.Contains("Default")) mode=0;
4092 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4094 Int_t shieldindex=0;
4095 if (shield.Contains("SDD")) { // SDDouter
4096 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4098 } else if (shield.Contains("SPD")) { // SPDouter
4099 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4102 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4106 // do nothing if we are already beyond the shield
4107 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4108 if(dir<0 && rTrack > rToGo) return 1; // going outward
4109 if(dir>0 && rTrack < rToGo) return 1; // going inward
4113 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4115 Int_t index=2*fCurrentEsdTrack+shieldindex;
4117 Double_t xOverX0,x0,lengthTimesMeanDensity;
4122 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4123 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4124 lengthTimesMeanDensity = xOverX0*x0;
4125 lengthTimesMeanDensity *= dir;
4126 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4129 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4130 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4133 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4134 xOverX0 = fxOverX0Shield[shieldindex];
4135 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4136 lengthTimesMeanDensity *= dir;
4137 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4140 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4141 if(fxOverX0ShieldTrks[index]<0) {
4142 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4143 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4144 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4145 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4146 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4149 xOverX0 = fxOverX0ShieldTrks[index];
4150 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4151 lengthTimesMeanDensity *= dir;
4152 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4158 //------------------------------------------------------------------------
4159 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4161 Double_t oldGlobXYZ[3],
4162 TString direction) {
4163 //-------------------------------------------------------------------
4164 // Propagate beyond layer and correct for material
4165 // (material budget in different ways according to fUseTGeo value)
4166 // Add time if going outward (PropagateTo or PropagateToTGeo)
4167 //-------------------------------------------------------------------
4169 // Define budget mode:
4170 // 0: material from AliITSRecoParam (hard coded)
4171 // 1: material from TGeo in stepsof X cm (on the fly)
4172 // X = AliITSRecoParam::GetStepSizeTGeo()
4173 // 2: material from lut
4174 // 3: material from TGeo in one step (same for all hypotheses)
4187 if(fTrackingPhase.Contains("Clusters2Tracks"))
4188 { mode=3; } else { mode=1; }
4191 if(fTrackingPhase.Contains("Clusters2Tracks"))
4192 { mode=3; } else { mode=2; }
4198 if(fTrackingPhase.Contains("Default")) mode=0;
4200 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4202 Double_t r=fgLayers[layerindex].GetR();
4203 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4205 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4207 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4209 Int_t index=6*fCurrentEsdTrack+layerindex;
4212 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4215 // back before material (no correction)
4217 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4218 if (!t->GetLocalXat(rOld,xOld)) return 0;
4219 if (!t->Propagate(xOld)) return 0;
4223 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4224 lengthTimesMeanDensity = xOverX0*x0;
4225 lengthTimesMeanDensity *= dir;
4226 // Bring the track beyond the material
4227 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4230 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4231 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4234 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4235 xOverX0 = fxOverX0Layer[layerindex];
4236 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4237 lengthTimesMeanDensity *= dir;
4238 // Bring the track beyond the material
4239 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4242 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4243 if(fxOverX0LayerTrks[index]<0) {
4244 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4245 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4246 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4247 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4248 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4249 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4252 xOverX0 = fxOverX0LayerTrks[index];
4253 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4254 lengthTimesMeanDensity *= dir;
4255 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4262 //------------------------------------------------------------------------
4263 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4264 //-----------------------------------------------------------------
4265 // Initialize LUT for storing material for each prolonged track
4266 //-----------------------------------------------------------------
4267 fxOverX0PipeTrks = new Float_t[ntracks];
4268 fxTimesRhoPipeTrks = new Float_t[ntracks];
4269 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4270 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4271 fxOverX0LayerTrks = new Float_t[ntracks*6];
4272 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4274 for(Int_t i=0; i<ntracks; i++) {
4275 fxOverX0PipeTrks[i] = -1.;
4276 fxTimesRhoPipeTrks[i] = -1.;
4278 for(Int_t j=0; j<ntracks*2; j++) {
4279 fxOverX0ShieldTrks[j] = -1.;
4280 fxTimesRhoShieldTrks[j] = -1.;
4282 for(Int_t k=0; k<ntracks*6; k++) {
4283 fxOverX0LayerTrks[k] = -1.;
4284 fxTimesRhoLayerTrks[k] = -1.;
4291 //------------------------------------------------------------------------
4292 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4293 //-----------------------------------------------------------------
4294 // Delete LUT for storing material for each prolonged track
4295 //-----------------------------------------------------------------
4296 if(fxOverX0PipeTrks) {
4297 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4299 if(fxOverX0ShieldTrks) {
4300 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4303 if(fxOverX0LayerTrks) {
4304 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4306 if(fxTimesRhoPipeTrks) {
4307 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4309 if(fxTimesRhoShieldTrks) {
4310 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4312 if(fxTimesRhoLayerTrks) {
4313 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4317 //------------------------------------------------------------------------
4318 void AliITStrackerMI::SetForceSkippingOfLayer() {
4319 //-----------------------------------------------------------------
4320 // Check if we are forced to skip layers
4321 // either we set to skip them in RecoParam
4322 // or they were off during data-taking
4323 //-----------------------------------------------------------------
4325 const AliEventInfo *eventInfo = GetEventInfo();
4327 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4328 fForceSkippingOfLayer[l] = 0;
4330 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4334 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4335 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4337 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4338 } else if(l==2 || l==3) {
4339 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4341 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4347 //------------------------------------------------------------------------
4348 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4349 Int_t ilayer,Int_t idet) const {
4350 //-----------------------------------------------------------------
4351 // This method is used to decide whether to allow a prolongation
4352 // without clusters, because we want to skip the layer.
4353 // In this case the return value is > 0:
4354 // return 1: the user requested to skip a layer
4355 // return 2: track outside z acceptance
4356 //-----------------------------------------------------------------
4358 if (ForceSkippingOfLayer(ilayer)) return 1;
4360 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4362 if (idet<0 && // out in z
4363 ilayer>innerLayCanSkip &&
4364 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4365 // check if track will cross SPD outer layer
4366 Double_t phiAtSPD2,zAtSPD2;
4367 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4368 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4370 return 2; // always allow skipping, changed on 05.11.2009
4375 //------------------------------------------------------------------------
4376 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4377 Int_t ilayer,Int_t idet,
4378 Double_t dz,Double_t dy,
4379 Bool_t noClusters) const {
4380 //-----------------------------------------------------------------
4381 // This method is used to decide whether to allow a prolongation
4382 // without clusters, because there is a dead zone in the road.
4383 // In this case the return value is > 0:
4384 // return 1: dead zone at z=0,+-7cm in SPD
4385 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4386 // return 2: all road is "bad" (dead or noisy) from the OCDB
4387 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4388 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4389 //-----------------------------------------------------------------
4391 // check dead zones at z=0,+-7cm in the SPD
4392 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4393 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4394 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4395 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4396 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4397 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4398 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4399 for (Int_t i=0; i<3; i++)
4400 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4401 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4402 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4406 // check bad zones from OCDB
4407 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4409 if (idet<0) return 0;
4411 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4414 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4415 if (ilayer==0 || ilayer==1) { // ---------- SPD
4417 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4419 detSizeFactorX *= 2.;
4420 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4423 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4424 if (detType==2) segm->SetLayer(ilayer+1);
4425 Float_t detSizeX = detSizeFactorX*segm->Dx();
4426 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4428 // check if the road overlaps with bad chips
4430 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
4431 Float_t zlocmin = zloc-dz;
4432 Float_t zlocmax = zloc+dz;
4433 Float_t xlocmin = xloc-dy;
4434 Float_t xlocmax = xloc+dy;
4435 Int_t chipsInRoad[100];
4437 // check if road goes out of detector
4438 Bool_t touchNeighbourDet=kFALSE;
4439 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4440 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4441 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4442 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4443 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));
4445 // check if this detector is bad
4447 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4448 if(!touchNeighbourDet) {
4449 return 2; // all detectors in road are bad
4451 return 3; // at least one is bad
4455 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4456 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4457 if (!nChipsInRoad) return 0;
4459 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4460 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4461 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4462 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4463 if (det.IsChipBad(chipsInRoad[iCh])) {
4471 if(!touchNeighbourDet) {
4472 AliDebug(2,"all bad in road");
4473 return 2; // all chips in road are bad
4475 return 3; // at least a bad chip in road
4480 AliDebug(2,"at least a bad in road");
4481 return 3; // at least a bad chip in road
4485 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4486 || !noClusters) return 0;
4488 // There are no clusters in road: check if there is at least
4489 // a bad SPD pixel or SDD anode or SSD strips on both sides
4491 Int_t idetInITS=idet;
4492 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4494 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4495 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4498 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4502 //------------------------------------------------------------------------
4503 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4504 const AliITStrackMI *track,
4505 Float_t &xloc,Float_t &zloc) const {
4506 //-----------------------------------------------------------------
4507 // Gives position of track in local module ref. frame
4508 //-----------------------------------------------------------------
4513 if(idet<0) return kTRUE; // track out of z acceptance of layer
4515 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4517 Int_t lad = Int_t(idet/ndet) + 1;
4519 Int_t det = idet - (lad-1)*ndet + 1;
4521 Double_t xyzGlob[3],xyzLoc[3];
4523 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4524 // take into account the misalignment: xyz at real detector plane
4525 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4527 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4529 xloc = (Float_t)xyzLoc[0];
4530 zloc = (Float_t)xyzLoc[2];
4534 //------------------------------------------------------------------------
4535 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4537 // Method to be optimized further:
4538 // Aim: decide whether a track can be used for PlaneEff evaluation
4539 // the decision is taken based on the track quality at the layer under study
4540 // no information on the clusters on this layer has to be used
4541 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4542 // the cut is done on number of sigmas from the boundaries
4544 // Input: Actual track, layer [0,5] under study
4546 // Return: kTRUE if this is a good track
4548 // it will apply a pre-selection to obtain good quality tracks.
4549 // Here also you will have the possibility to put a control on the
4550 // impact point of the track on the basic block, in order to exclude border regions
4551 // this will be done by calling a proper method of the AliITSPlaneEff class.
4553 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4554 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4556 Int_t index[AliITSgeomTGeo::kNLayers];
4558 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4560 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4561 index[k]=clusters[k];
4565 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4566 AliITSlayer &layer=fgLayers[ilayer];
4567 Double_t r=layer.GetR();
4568 AliITStrackMI tmp(*track);
4570 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4572 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
4573 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4574 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4575 if (tmp.GetClIndex(lay)>=0) ncl++;
4577 Bool_t nextout = kFALSE;
4578 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4579 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4580 Bool_t nextin = kFALSE;
4581 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4582 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4583 if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4585 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4586 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4587 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4588 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4592 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4593 Int_t idet=layer.FindDetectorIndex(phi,z);
4594 if(idet<0) { AliInfo(Form("cannot find detector"));
4597 // here check if it has good Chi Square.
4599 //propagate to the intersection with the detector plane
4600 const AliITSdetector &det=layer.GetDetector(idet);
4601 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4605 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4606 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4607 if(key>fPlaneEff->Nblock()) return kFALSE;
4608 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4609 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4611 // DEFINITION OF SEARCH ROAD FOR accepting a track
4613 //For the time being they are hard-wired, later on from AliITSRecoParam
4614 // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
4615 // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
4618 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4619 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4620 // done for RecPoints
4622 // exclude tracks at boundary between detectors
4623 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4624 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4625 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4626 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4627 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4629 if ( (locx-dx < blockXmn+boundaryWidth) ||
4630 (locx+dx > blockXmx-boundaryWidth) ||
4631 (locz-dz < blockZmn+boundaryWidth) ||
4632 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4635 //------------------------------------------------------------------------
4636 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4638 // This Method has to be optimized! For the time-being it uses the same criteria
4639 // as those used in the search of extra clusters for overlapping modules.
4641 // Method Purpose: estabilish whether a track has produced a recpoint or not
4642 // in the layer under study (For Plane efficiency)
4644 // inputs: AliITStrackMI* track (pointer to a usable track)
4646 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4647 // traversed by this very track. In details:
4648 // - if a cluster can be associated to the track then call
4649 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4650 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4653 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4654 AliITSlayer &layer=fgLayers[ilayer];
4655 Double_t r=layer.GetR();
4656 AliITStrackMI tmp(*track);
4660 if (!tmp.GetPhiZat(r,phi,z)) return;
4661 Int_t idet=layer.FindDetectorIndex(phi,z);
4663 if(idet<0) { AliInfo(Form("cannot find detector"));
4667 //propagate to the intersection with the detector plane
4668 const AliITSdetector &det=layer.GetDetector(idet);
4669 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4673 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4675 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4676 TMath::Sqrt(tmp.GetSigmaZ2() +
4677 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4678 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4679 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4680 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4681 TMath::Sqrt(tmp.GetSigmaY2() +
4682 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4683 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4684 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4686 // road in global (rphi,z) [i.e. in tracking ref. system]
4687 Double_t zmin = tmp.GetZ() - dz;
4688 Double_t zmax = tmp.GetZ() + dz;
4689 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4690 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4692 // select clusters in road
4693 layer.SelectClusters(zmin,zmax,ymin,ymax);
4695 // Define criteria for track-cluster association
4696 Double_t msz = tmp.GetSigmaZ2() +
4697 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4698 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4699 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4700 Double_t msy = tmp.GetSigmaY2() +
4701 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4702 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4703 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4704 if (tmp.GetConstrain()) {
4705 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4706 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4708 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4709 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4711 msz = 1./msz; // 1/RoadZ^2
4712 msy = 1./msy; // 1/RoadY^2
4715 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4717 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4718 //Double_t tolerance=0.2;
4719 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4720 idetc = cl->GetDetectorIndex();
4721 if(idet!=idetc) continue;
4722 //Int_t ilay = cl->GetLayer();
4724 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4725 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4727 Double_t chi2=tmp.GetPredictedChi2(cl);
4728 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4732 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4734 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4735 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4736 if(key>fPlaneEff->Nblock()) return;
4737 Bool_t found=kFALSE;
4740 while ((cl=layer.GetNextCluster(clidx))!=0) {
4741 idetc = cl->GetDetectorIndex();
4742 if(idet!=idetc) continue;
4743 // here real control to see whether the cluster can be associated to the track.
4744 // cluster not associated to track
4745 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4746 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
4747 // calculate track-clusters chi2
4748 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4749 // in particular, the error associated to the cluster
4750 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4752 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4754 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4755 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4756 // track->SetExtraModule(ilayer,idetExtra);
4758 if(!fPlaneEff->UpDatePlaneEff(found,key))
4759 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4760 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4761 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
4762 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4763 Int_t cltype[2]={-999,-999};
4767 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4768 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4771 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4772 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4773 cltype[0]=layer.GetCluster(ci)->GetNy();
4774 cltype[1]=layer.GetCluster(ci)->GetNz();
4776 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4777 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
4778 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4779 // It is computed properly by calling the method
4780 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4782 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4783 //if (tmp.PropagateTo(x,0.,0.)) {
4784 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4785 AliCluster c(*layer.GetCluster(ci));
4786 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4787 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4788 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
4789 clu[2]=TMath::Sqrt(c.GetSigmaY2());
4790 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4793 fPlaneEff->FillHistos(key,found,tr,clu,cltype);