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>
38 #include "AliGeomManager.h"
39 #include "AliITSPlaneEff.h"
40 #include "AliTrackPointArray.h"
41 #include "AliESDEvent.h"
42 #include "AliESDtrack.h"
44 #include "AliITSChannelStatus.h"
45 #include "AliITSDetTypeRec.h"
46 #include "AliITSRecPoint.h"
47 #include "AliITSRecPointContainer.h"
48 #include "AliITSgeomTGeo.h"
49 #include "AliITSReconstructor.h"
50 #include "AliITSClusterParam.h"
51 #include "AliITSsegmentation.h"
52 #include "AliITSCalibration.h"
53 #include "AliITSPlaneEffSPD.h"
54 #include "AliITSPlaneEffSDD.h"
55 #include "AliITSPlaneEffSSD.h"
56 #include "AliITSV0Finder.h"
57 #include "AliITStrackerMI.h"
58 #include "AliMathBase.h"
60 ClassImp(AliITStrackerMI)
62 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
64 AliITStrackerMI::AliITStrackerMI():AliTracker(),
74 fLastLayerToTrackTo(0),
77 fTrackingPhase("Default"),
83 fxTimesRhoPipeTrks(0),
84 fxOverX0ShieldTrks(0),
85 fxTimesRhoShieldTrks(0),
87 fxTimesRhoLayerTrks(0),
94 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
96 fxOverX0Shield[i]=-1.;
97 fxTimesRhoShield[i]=-1.;
100 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
101 fOriginal.SetOwner();
102 for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
103 for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
106 //------------------------------------------------------------------------
107 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
108 fI(AliITSgeomTGeo::GetNLayers()),
117 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
120 fTrackingPhase("Default"),
126 fxTimesRhoPipeTrks(0),
127 fxOverX0ShieldTrks(0),
128 fxTimesRhoShieldTrks(0),
129 fxOverX0LayerTrks(0),
130 fxTimesRhoLayerTrks(0),
132 fITSChannelStatus(0),
135 //--------------------------------------------------------------------
136 //This is the AliITStrackerMI constructor
137 //--------------------------------------------------------------------
139 AliWarning("\"geom\" is actually a dummy argument !");
142 fOriginal.SetOwner();
146 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
147 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
148 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
150 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
151 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
152 Double_t poff=TMath::ATan2(y,x);
155 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
156 Double_t r=TMath::Sqrt(x*x + y*y);
158 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
159 r += TMath::Sqrt(x*x + y*y);
160 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
161 r += TMath::Sqrt(x*x + y*y);
162 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
163 r += TMath::Sqrt(x*x + y*y);
166 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
168 for (Int_t j=1; j<nlad+1; j++) {
169 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
170 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
171 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
173 Double_t txyz[3]={0.};
174 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
175 m.LocalToMaster(txyz,xyz);
176 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
177 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
179 if (phi<0) phi+=TMath::TwoPi();
180 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
182 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
183 new(&det) AliITSdetector(r,phi);
184 // compute the real radius (with misalignment)
185 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
187 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
188 mmisal.LocalToMaster(txyz,xyz);
189 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
190 det.SetRmisal(rmisal);
192 } // end loop on detectors
193 } // end loop on ladders
194 fForceSkippingOfLayer[i-1] = 0;
195 } // end loop on layers
198 fI=AliITSgeomTGeo::GetNLayers();
201 fConstraint[0]=1; fConstraint[1]=0;
203 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
204 AliITSReconstructor::GetRecoParam()->GetYVdef(),
205 AliITSReconstructor::GetRecoParam()->GetZVdef()};
206 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
207 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
208 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
209 SetVertex(xyzVtx,ersVtx);
211 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
212 for (Int_t i=0;i<100000;i++){
213 fBestTrackIndex[i]=0;
216 // store positions of centre of SPD modules (in z)
217 // The convetion is that fSPDdetzcentre is ordered from -z to +z
219 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
220 if (tr[2]<0) { // old geom (up to v5asymmPPR)
221 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
222 fSPDdetzcentre[0] = tr[2];
223 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
224 fSPDdetzcentre[1] = tr[2];
225 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
226 fSPDdetzcentre[2] = tr[2];
227 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
228 fSPDdetzcentre[3] = tr[2];
229 } else { // new geom (from v11Hybrid)
230 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
231 fSPDdetzcentre[0] = tr[2];
232 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
233 fSPDdetzcentre[1] = tr[2];
234 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
235 fSPDdetzcentre[2] = tr[2];
236 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
237 fSPDdetzcentre[3] = tr[2];
240 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
241 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
242 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
246 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
247 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
249 if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
250 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
252 // only for plane efficiency evaluation
253 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
254 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
255 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
256 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
257 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
258 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
259 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
260 else fPlaneEff = new AliITSPlaneEffSSD();
261 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
262 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
263 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
267 //------------------------------------------------------------------------
268 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
270 fBestTrack(tracker.fBestTrack),
271 fTrackToFollow(tracker.fTrackToFollow),
272 fTrackHypothesys(tracker.fTrackHypothesys),
273 fBestHypothesys(tracker.fBestHypothesys),
274 fOriginal(tracker.fOriginal),
275 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
276 fPass(tracker.fPass),
277 fAfterV0(tracker.fAfterV0),
278 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
279 fCoefficients(tracker.fCoefficients),
281 fTrackingPhase(tracker.fTrackingPhase),
282 fUseTGeo(tracker.fUseTGeo),
283 fNtracks(tracker.fNtracks),
284 fxOverX0Pipe(tracker.fxOverX0Pipe),
285 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
287 fxTimesRhoPipeTrks(0),
288 fxOverX0ShieldTrks(0),
289 fxTimesRhoShieldTrks(0),
290 fxOverX0LayerTrks(0),
291 fxTimesRhoLayerTrks(0),
292 fDebugStreamer(tracker.fDebugStreamer),
293 fITSChannelStatus(tracker.fITSChannelStatus),
294 fkDetTypeRec(tracker.fkDetTypeRec),
295 fPlaneEff(tracker.fPlaneEff) {
297 fOriginal.SetOwner();
300 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
303 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
304 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
307 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
308 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
312 //------------------------------------------------------------------------
313 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
314 //Assignment operator
315 this->~AliITStrackerMI();
316 new(this) AliITStrackerMI(tracker);
320 //------------------------------------------------------------------------
321 AliITStrackerMI::~AliITStrackerMI()
326 if (fCoefficients) delete [] fCoefficients;
327 DeleteTrksMaterialLUT();
328 if (fDebugStreamer) {
329 //fDebugStreamer->Close();
330 delete fDebugStreamer;
332 if(fITSChannelStatus) delete fITSChannelStatus;
333 if(fPlaneEff) delete fPlaneEff;
335 //------------------------------------------------------------------------
336 void AliITStrackerMI::ReadBadFromDetTypeRec() {
337 //--------------------------------------------------------------------
338 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
340 //--------------------------------------------------------------------
342 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
344 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
346 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
349 if(fITSChannelStatus) delete fITSChannelStatus;
350 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
352 // ITS detectors and chips
353 Int_t i=0,j=0,k=0,ndet=0;
354 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
355 Int_t nBadDetsPerLayer=0;
356 ndet=AliITSgeomTGeo::GetNDetectors(i);
357 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
358 for (k=1; k<ndet+1; k++) {
359 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
360 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
361 if(det.IsBad()) {nBadDetsPerLayer++;}
362 } // end loop on detectors
363 } // end loop on ladders
364 AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
365 } // end loop on layers
369 //------------------------------------------------------------------------
370 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
371 //--------------------------------------------------------------------
372 //This function loads ITS clusters
373 //--------------------------------------------------------------------
375 TClonesArray *clusters = NULL;
376 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
377 clusters=rpcont->FetchClusters(0,cTree);
378 if(!clusters) return 1;
380 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
381 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
384 Int_t i=0,j=0,ndet=0;
386 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
387 ndet=fgLayers[i].GetNdetectors();
388 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
389 for (; j<jmax; j++) {
390 // if (!cTree->GetEvent(j)) continue;
391 clusters = rpcont->UncheckedGetClusters(j);
392 if(!clusters)continue;
393 Int_t ncl=clusters->GetEntriesFast();
394 SignDeltas(clusters,GetZ());
397 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
398 detector=c->GetDetectorIndex();
400 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
402 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
404 AliWarning(Form("Too many clusters on layer %d!",i));
409 // add dead zone "virtual" cluster in SPD, if there is a cluster within
410 // zwindow cm from the dead zone
411 // This method assumes that fSPDdetzcentre is ordered from -z to +z
412 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
413 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
414 Int_t lab[4] = {0,0,0,detector};
415 Int_t info[3] = {0,0,i};
416 Float_t q = 0.; // this identifies virtual clusters
417 Float_t hit[6] = {xdead,
419 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
420 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
423 Bool_t local = kTRUE;
424 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
425 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
426 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
427 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
428 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
429 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
430 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
431 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
432 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
433 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
434 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
435 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
436 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
437 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
438 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
439 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
440 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
441 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
442 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
444 } // "virtual" clusters in SPD
448 fgLayers[i].ResetRoad(); //road defined by the cluster density
449 fgLayers[i].SortClusters();
452 // check whether we have to skip some layers
453 SetForceSkippingOfLayer();
457 //------------------------------------------------------------------------
458 void AliITStrackerMI::UnloadClusters() {
459 //--------------------------------------------------------------------
460 //This function unloads ITS clusters
461 //--------------------------------------------------------------------
462 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
464 //------------------------------------------------------------------------
465 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
466 //--------------------------------------------------------------------
467 // Publishes all pointers to clusters known to the tracker into the
468 // passed object array.
469 // The ownership is not transfered - the caller is not expected to delete
471 //--------------------------------------------------------------------
473 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
474 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
475 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
482 //------------------------------------------------------------------------
483 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
484 //--------------------------------------------------------------------
485 // Correction for the material between the TPC and the ITS
486 //--------------------------------------------------------------------
487 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
488 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
489 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
490 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
491 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
492 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
493 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
494 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
496 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
502 //------------------------------------------------------------------------
503 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
504 //--------------------------------------------------------------------
505 // This functions reconstructs ITS tracks
506 // The clusters must be already loaded !
507 //--------------------------------------------------------------------
509 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
511 fTrackingPhase="Clusters2Tracks";
513 TObjArray itsTracks(15000);
515 fEsd = event; // store pointer to the esd
517 // temporary (for cosmics)
518 if(event->GetVertex()) {
519 TString title = event->GetVertex()->GetTitle();
520 if(title.Contains("cosmics")) {
521 Double_t xyz[3]={GetX(),GetY(),GetZ()};
522 Double_t exyz[3]={0.1,0.1,0.1};
528 {/* Read ESD tracks */
529 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
530 Int_t nentr=event->GetNumberOfTracks();
532 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
534 AliESDtrack *esd=event->GetTrack(nentr);
535 // ---- for debugging:
536 //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); }
538 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
539 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
540 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
541 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
542 AliITStrackMI *t = new AliITStrackMI(*esd);
543 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
544 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
547 // look at the ESD mass hypothesys !
548 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
550 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
552 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
553 //track - can be V0 according to TPC
555 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
559 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
563 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
567 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
572 t->SetReconstructed(kFALSE);
573 itsTracks.AddLast(t);
574 fOriginal.AddLast(t);
576 } /* End Read ESD tracks */
580 Int_t nentr=itsTracks.GetEntriesFast();
581 fTrackHypothesys.Expand(nentr);
582 fBestHypothesys.Expand(nentr);
583 MakeCoefficients(nentr);
584 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
586 // THE TWO TRACKING PASSES
587 for (fPass=0; fPass<2; fPass++) {
588 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
589 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
590 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
591 if (t==0) continue; //this track has been already tracked
592 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
593 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
594 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
595 if (fConstraint[fPass]) {
596 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
597 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
600 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
601 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
603 ResetTrackToFollow(*t);
606 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
609 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
611 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
612 if (!besttrack) continue;
613 besttrack->SetLabel(tpcLabel);
614 // besttrack->CookdEdx();
616 besttrack->SetFakeRatio(1.);
617 CookLabel(besttrack,0.); //For comparison only
618 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
620 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
622 t->SetReconstructed(kTRUE);
624 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
626 GetBestHypothesysMIP(itsTracks);
627 } // end loop on the two tracking passes
629 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
630 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
635 Int_t entries = fTrackHypothesys.GetEntriesFast();
636 for (Int_t ientry=0; ientry<entries; ientry++) {
637 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
638 if (array) array->Delete();
639 delete fTrackHypothesys.RemoveAt(ientry);
642 fTrackHypothesys.Delete();
643 entries = fBestHypothesys.GetEntriesFast();
644 for (Int_t ientry=0; ientry<entries; ientry++) {
645 TObjArray * array =(TObjArray*)fBestHypothesys.UncheckedAt(ientry);
646 if (array) array->Delete();
647 delete fBestHypothesys.RemoveAt(ientry);
649 fBestHypothesys.Delete();
651 delete [] fCoefficients;
653 DeleteTrksMaterialLUT();
655 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
657 fTrackingPhase="Default";
661 //------------------------------------------------------------------------
662 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
663 //--------------------------------------------------------------------
664 // This functions propagates reconstructed ITS tracks back
665 // The clusters must be loaded !
666 //--------------------------------------------------------------------
667 fTrackingPhase="PropagateBack";
668 Int_t nentr=event->GetNumberOfTracks();
669 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
672 for (Int_t i=0; i<nentr; i++) {
673 AliESDtrack *esd=event->GetTrack(i);
675 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
676 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
678 AliITStrackMI *t = new AliITStrackMI(*esd);
680 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
682 ResetTrackToFollow(*t);
685 // propagate to vertex [SR, GSI 17.02.2003]
686 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
687 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
688 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
689 fTrackToFollow.StartTimeIntegral();
690 // from vertex to outside pipe
691 CorrectForPipeMaterial(&fTrackToFollow,"outward");
693 // Start time integral and add distance from current position to vertex
694 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
695 fTrackToFollow.GetXYZ(xyzTrk);
697 for (Int_t icoord=0; icoord<3; icoord++) {
698 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
701 fTrackToFollow.StartTimeIntegral();
702 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
704 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
705 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
706 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
710 fTrackToFollow.SetLabel(t->GetLabel());
711 //fTrackToFollow.CookdEdx();
712 CookLabel(&fTrackToFollow,0.); //For comparison only
713 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
714 //UseClusters(&fTrackToFollow);
720 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
722 fTrackingPhase="Default";
726 //------------------------------------------------------------------------
727 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
728 //--------------------------------------------------------------------
729 // This functions refits ITS tracks using the
730 // "inward propagated" TPC tracks
731 // The clusters must be loaded !
732 //--------------------------------------------------------------------
733 fTrackingPhase="RefitInward";
735 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
737 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
738 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
740 Int_t nentr=event->GetNumberOfTracks();
741 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
744 for (Int_t i=0; i<nentr; i++) {
745 AliESDtrack *esd=event->GetTrack(i);
747 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
748 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
749 if (esd->GetStatus()&AliESDtrack::kTPCout)
750 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
752 AliITStrackMI *t = new AliITStrackMI(*esd);
754 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
755 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
760 ResetTrackToFollow(*t);
761 fTrackToFollow.ResetClusters();
763 // ITS standalone tracks
764 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
765 fTrackToFollow.ResetCovariance(10.);
766 // protection for loopers that can have parameters screwed up
767 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
768 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
775 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
776 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
778 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
779 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
780 AliDebug(2," refit OK");
781 fTrackToFollow.SetLabel(t->GetLabel());
782 // fTrackToFollow.CookdEdx();
783 CookdEdx(&fTrackToFollow);
785 CookLabel(&fTrackToFollow,0.0); //For comparison only
788 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
789 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
790 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
791 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
792 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
793 Double_t r[3]={0.,0.,0.};
795 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
802 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
804 fTrackingPhase="Default";
808 //------------------------------------------------------------------------
809 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
810 //--------------------------------------------------------------------
811 // Return pointer to a given cluster
812 //--------------------------------------------------------------------
813 Int_t l=(index & 0xf0000000) >> 28;
814 Int_t c=(index & 0x0fffffff) >> 00;
815 return fgLayers[l].GetCluster(c);
817 //------------------------------------------------------------------------
818 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
819 //--------------------------------------------------------------------
820 // Get track space point with index i
821 //--------------------------------------------------------------------
823 Int_t l=(index & 0xf0000000) >> 28;
824 Int_t c=(index & 0x0fffffff) >> 00;
825 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
826 Int_t idet = cl->GetDetectorIndex();
830 cl->GetGlobalXYZ(xyz);
831 cl->GetGlobalCov(cov);
833 p.SetCharge(cl->GetQ());
834 p.SetDriftTime(cl->GetDriftTime());
835 p.SetChargeRatio(cl->GetChargeRatio());
836 p.SetClusterType(cl->GetClusterType());
837 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
840 iLayer = AliGeomManager::kSPD1;
843 iLayer = AliGeomManager::kSPD2;
846 iLayer = AliGeomManager::kSDD1;
849 iLayer = AliGeomManager::kSDD2;
852 iLayer = AliGeomManager::kSSD1;
855 iLayer = AliGeomManager::kSSD2;
858 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
861 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
862 p.SetVolumeID((UShort_t)volid);
865 //------------------------------------------------------------------------
866 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
867 AliTrackPoint& p, const AliESDtrack *t) {
868 //--------------------------------------------------------------------
869 // Get track space point with index i
870 // (assign error estimated during the tracking)
871 //--------------------------------------------------------------------
873 Int_t l=(index & 0xf0000000) >> 28;
874 Int_t c=(index & 0x0fffffff) >> 00;
875 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
876 Int_t idet = cl->GetDetectorIndex();
878 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
880 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
882 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
883 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
884 Double_t alpha = t->GetAlpha();
885 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
886 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
887 phi += alpha-det.GetPhi();
888 Float_t tgphi = TMath::Tan(phi);
890 Float_t tgl = t->GetTgl(); // tgl about const along track
891 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
893 Float_t errtrky,errtrkz,covyz;
894 Bool_t addMisalErr=kFALSE;
895 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
899 cl->GetGlobalXYZ(xyz);
900 // cl->GetGlobalCov(cov);
901 Float_t pos[3] = {0.,0.,0.};
902 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
903 tmpcl.GetGlobalCov(cov);
906 p.SetCharge(cl->GetQ());
907 p.SetDriftTime(cl->GetDriftTime());
908 p.SetChargeRatio(cl->GetChargeRatio());
909 p.SetClusterType(cl->GetClusterType());
911 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
914 iLayer = AliGeomManager::kSPD1;
917 iLayer = AliGeomManager::kSPD2;
920 iLayer = AliGeomManager::kSDD1;
923 iLayer = AliGeomManager::kSDD2;
926 iLayer = AliGeomManager::kSSD1;
929 iLayer = AliGeomManager::kSSD2;
932 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
935 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
937 p.SetVolumeID((UShort_t)volid);
940 //------------------------------------------------------------------------
941 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
943 //--------------------------------------------------------------------
944 // Follow prolongation tree
945 //--------------------------------------------------------------------
947 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
948 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
951 AliESDtrack * esd = otrack->GetESDtrack();
952 if (esd->GetV0Index(0)>0) {
953 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
954 // mapping of ESD track is different as ITS track in Containers
955 // Need something more stable
956 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
957 for (Int_t i=0;i<3;i++){
958 Int_t index = esd->GetV0Index(i);
960 AliESDv0 * vertex = fEsd->GetV0(index);
961 if (vertex->GetStatus()<0) continue; // rejected V0
963 if (esd->GetSign()>0) {
964 vertex->SetIndex(0,esdindex);
966 vertex->SetIndex(1,esdindex);
970 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
972 bestarray = new TObjArray(5);
973 bestarray->SetOwner();
974 fBestHypothesys.AddAt(bestarray,esdindex);
978 //setup tree of the prolongations
980 static AliITStrackMI tracks[7][100];
981 AliITStrackMI *currenttrack;
982 static AliITStrackMI currenttrack1;
983 static AliITStrackMI currenttrack2;
984 static AliITStrackMI backuptrack;
986 Int_t nindexes[7][100];
987 Float_t normalizedchi2[100];
988 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
989 otrack->SetNSkipped(0);
990 new (&(tracks[6][0])) AliITStrackMI(*otrack);
992 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
993 Int_t modstatus = 1; // found
997 // follow prolongations
998 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
999 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1002 AliITSlayer &layer=fgLayers[ilayer];
1003 Double_t r = layer.GetR();
1009 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1011 if (ntracks[ilayer]>=100) break;
1012 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1013 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1014 if (ntracks[ilayer]>15+ilayer){
1015 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1016 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1019 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1021 // material between SSD and SDD, SDD and SPD
1023 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1025 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1029 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1030 Int_t idet=layer.FindDetectorIndex(phi,z);
1032 Double_t trackGlobXYZ1[3];
1033 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1035 // Get the budget to the primary vertex for the current track being prolonged
1036 Double_t budgetToPrimVertex = GetEffectiveThickness();
1038 // check if we allow a prolongation without point
1039 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1041 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1042 // propagate to the layer radius
1043 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1044 if(!vtrack->Propagate(xToGo)) continue;
1045 // apply correction for material of the current layer
1046 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1047 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1048 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1049 vtrack->SetClIndex(ilayer,-1);
1050 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1051 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1052 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1054 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1059 // track outside layer acceptance in z
1060 if (idet<0) continue;
1062 //propagate to the intersection with the detector plane
1063 const AliITSdetector &det=layer.GetDetector(idet);
1064 new(¤ttrack2) AliITStrackMI(currenttrack1);
1065 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1066 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1067 currenttrack1.SetDetectorIndex(idet);
1068 currenttrack2.SetDetectorIndex(idet);
1069 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1072 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1074 // road in global (rphi,z) [i.e. in tracking ref. system]
1075 Double_t zmin,zmax,ymin,ymax;
1076 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1078 // select clusters in road
1079 layer.SelectClusters(zmin,zmax,ymin,ymax);
1080 //********************
1082 // Define criteria for track-cluster association
1083 Double_t msz = currenttrack1.GetSigmaZ2() +
1084 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1085 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1086 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1087 Double_t msy = currenttrack1.GetSigmaY2() +
1088 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1089 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1090 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1092 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1093 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1095 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1096 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1098 msz = 1./msz; // 1/RoadZ^2
1099 msy = 1./msy; // 1/RoadY^2
1103 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1105 const AliITSRecPoint *cl=0;
1107 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1108 Bool_t deadzoneSPD=kFALSE;
1109 currenttrack = ¤ttrack1;
1111 // check if the road contains a dead zone
1112 Bool_t noClusters = kFALSE;
1113 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1114 if (noClusters) AliDebug(2,"no clusters in road");
1115 Double_t dz=0.5*(zmax-zmin);
1116 Double_t dy=0.5*(ymax-ymin);
1117 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1118 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1119 // create a prolongation without clusters (check also if there are no clusters in the road)
1122 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1123 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1124 updatetrack->SetClIndex(ilayer,-1);
1126 modstatus = 5; // no cls in road
1127 } else if (dead==1) {
1128 modstatus = 7; // holes in z in SPD
1129 } else if (dead==2 || dead==3 || dead==4) {
1130 modstatus = 2; // dead from OCDB
1132 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1133 // apply correction for material of the current layer
1134 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1135 if (constrain) { // apply vertex constrain
1136 updatetrack->SetConstrain(constrain);
1137 Bool_t isPrim = kTRUE;
1138 if (ilayer<4) { // check that it's close to the vertex
1139 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1140 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1141 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1142 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1143 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1145 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1147 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1149 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1150 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1152 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1153 updatetrack->SetDeadZoneProbability(ilayer,1.);
1154 } else if (dead==4) { // at least a single dead channel from OCDB
1155 updatetrack->SetDeadZoneProbability(ilayer,0.);
1162 // loop over clusters in the road
1163 while ((cl=layer.GetNextCluster(clidx))!=0) {
1164 if (ntracks[ilayer]>95) break; //space for skipped clusters
1165 Bool_t changedet =kFALSE;
1166 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1167 Int_t idetc=cl->GetDetectorIndex();
1169 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1170 // take into account misalignment (bring track to real detector plane)
1171 Double_t xTrOrig = currenttrack->GetX();
1172 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1173 // a first cut on track-cluster distance
1174 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1175 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1176 { // cluster not associated to track
1177 AliDebug(2,"not associated");
1178 // MvL: added here as well
1179 // bring track back to ideal detector plane
1180 currenttrack->Propagate(xTrOrig);
1183 // bring track back to ideal detector plane
1184 if (!currenttrack->Propagate(xTrOrig)) continue;
1185 } else { // have to move track to cluster's detector
1186 const AliITSdetector &detc=layer.GetDetector(idetc);
1187 // a first cut on track-cluster distance
1189 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1190 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1191 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1192 continue; // cluster not associated to track
1194 new (&backuptrack) AliITStrackMI(currenttrack2);
1196 currenttrack =¤ttrack2;
1197 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1198 new (currenttrack) AliITStrackMI(backuptrack);
1202 currenttrack->SetDetectorIndex(idetc);
1203 // Get again the budget to the primary vertex
1204 // for the current track being prolonged, if had to change detector
1205 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1208 // calculate track-clusters chi2
1209 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1211 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1212 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1213 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1214 if (ntracks[ilayer]>=100) continue;
1215 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1216 updatetrack->SetClIndex(ilayer,-1);
1217 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1219 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1220 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1221 AliDebug(2,"update failed");
1224 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1225 modstatus = 1; // found
1226 } else { // virtual cluster in dead zone
1227 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1228 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1229 modstatus = 7; // holes in z in SPD
1233 Float_t xlocnewdet,zlocnewdet;
1234 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1235 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1238 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1240 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1242 // apply correction for material of the current layer
1243 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1245 if (constrain) { // apply vertex constrain
1246 updatetrack->SetConstrain(constrain);
1247 Bool_t isPrim = kTRUE;
1248 if (ilayer<4) { // check that it's close to the vertex
1249 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1250 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1251 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1252 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1253 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1255 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1256 } //apply vertex constrain
1258 } // create new hypothesis
1260 AliDebug(2,"chi2 too large");
1263 } // loop over possible prolongations
1265 // allow one prolongation without clusters
1266 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1267 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1268 // apply correction for material of the current layer
1269 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1270 vtrack->SetClIndex(ilayer,-1);
1271 modstatus = 3; // skipped
1272 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1273 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1274 vtrack->IncrementNSkipped();
1279 } // loop over tracks in layer ilayer+1
1281 //loop over track candidates for the current layer
1287 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1288 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1289 if (normalizedchi2[itrack] <
1290 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1294 if (constrain) { // constrain
1295 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1297 } else { // !constrain
1298 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1303 // sort tracks by increasing normalized chi2
1304 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1305 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1306 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1307 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1308 } // end loop over layers
1312 // Now select tracks to be kept
1314 Int_t max = constrain ? 20 : 5;
1316 // tracks that reach layer 0 (SPD inner)
1317 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1318 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1319 if (track.GetNumberOfClusters()<2) continue;
1320 if (!constrain && track.GetNormChi2(0) >
1321 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1324 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1327 // tracks that reach layer 1 (SPD outer)
1328 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1329 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1330 if (track.GetNumberOfClusters()<4) continue;
1331 if (!constrain && track.GetNormChi2(1) >
1332 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1333 if (constrain) track.IncrementNSkipped();
1335 track.SetD(0,track.GetD(GetX(),GetY()));
1336 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1337 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1338 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1341 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1344 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1346 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1347 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1348 if (track.GetNumberOfClusters()<3) continue;
1349 if (track.GetNormChi2(2) >
1350 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1351 track.SetD(0,track.GetD(GetX(),GetY()));
1352 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1353 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1354 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1356 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1362 // register best track of each layer - important for V0 finder
1364 for (Int_t ilayer=0;ilayer<5;ilayer++){
1365 if (ntracks[ilayer]==0) continue;
1366 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1367 if (track.GetNumberOfClusters()<1) continue;
1368 CookLabel(&track,0);
1369 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1373 // update TPC V0 information
1375 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1376 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1377 for (Int_t i=0;i<3;i++){
1378 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1379 if (index==0) break;
1380 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1381 if (vertex->GetStatus()<0) continue; // rejected V0
1383 if (otrack->GetSign()>0) {
1384 vertex->SetIndex(0,esdindex);
1387 vertex->SetIndex(1,esdindex);
1389 //find nearest layer with track info
1390 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1391 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1392 Int_t nearest = nearestold;
1393 for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1394 if (ntracks[nearest]==0){
1399 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1400 if (nearestold<5&&nearest<5){
1401 Bool_t accept = track.GetNormChi2(nearest)<10;
1403 if (track.GetSign()>0) {
1404 vertex->SetParamP(track);
1405 vertex->Update(fprimvertex);
1406 //vertex->SetIndex(0,track.fESDtrack->GetID());
1407 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1409 vertex->SetParamN(track);
1410 vertex->Update(fprimvertex);
1411 //vertex->SetIndex(1,track.fESDtrack->GetID());
1412 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1414 vertex->SetStatus(vertex->GetStatus()+1);
1416 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1423 //------------------------------------------------------------------------
1424 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1426 //--------------------------------------------------------------------
1429 return fgLayers[layer];
1431 //------------------------------------------------------------------------
1432 AliITStrackerMI::AliITSlayer::AliITSlayer():
1462 //--------------------------------------------------------------------
1463 //default AliITSlayer constructor
1464 //--------------------------------------------------------------------
1465 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1466 fClusterWeight[i]=0;
1467 fClusterTracks[0][i]=-1;
1468 fClusterTracks[1][i]=-1;
1469 fClusterTracks[2][i]=-1;
1470 fClusterTracks[3][i]=-1;
1477 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
1478 for (Int_t j1=0; j1<6; j1++) {
1479 fClusters5[j1][j]=0;
1480 fClusterIndex5[j1][j]=-1;
1489 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
1490 for (Int_t j1=0; j1<11; j1++) {
1491 fClusters10[j1][j]=0;
1492 fClusterIndex10[j1][j]=-1;
1501 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
1502 for (Int_t j1=0; j1<21; j1++) {
1503 fClusters20[j1][j]=0;
1504 fClusterIndex20[j1][j]=-1;
1512 for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
1517 //------------------------------------------------------------------------
1518 AliITStrackerMI::AliITSlayer::
1519 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1548 //--------------------------------------------------------------------
1549 //main AliITSlayer constructor
1550 //--------------------------------------------------------------------
1551 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1552 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1554 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1555 fClusterWeight[i]=0;
1556 fClusterTracks[0][i]=-1;
1557 fClusterTracks[1][i]=-1;
1558 fClusterTracks[2][i]=-1;
1559 fClusterTracks[3][i]=-1;
1567 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer5; j++) {
1568 for (Int_t j1=0; j1<6; j1++) {
1569 fClusters5[j1][j]=0;
1570 fClusterIndex5[j1][j]=-1;
1579 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer10; j++) {
1580 for (Int_t j1=0; j1<11; j1++) {
1581 fClusters10[j1][j]=0;
1582 fClusterIndex10[j1][j]=-1;
1591 for (Int_t j=0; j<AliITSRecoParam::fgkMaxClusterPerLayer20; j++) {
1592 for (Int_t j1=0; j1<21; j1++) {
1593 fClusters20[j1][j]=0;
1594 fClusterIndex20[j1][j]=-1;
1602 for(Int_t i=0;i<AliITSRecoParam::fgkMaxClusterPerLayer;i++){
1608 //------------------------------------------------------------------------
1609 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1611 fPhiOffset(layer.fPhiOffset),
1612 fNladders(layer.fNladders),
1613 fZOffset(layer.fZOffset),
1614 fNdetectors(layer.fNdetectors),
1615 fDetectors(layer.fDetectors),
1620 fClustersCs(layer.fClustersCs),
1621 fClusterIndexCs(layer.fClusterIndexCs),
1625 fCurrentSlice(layer.fCurrentSlice),
1633 fAccepted(layer.fAccepted),
1635 fMaxSigmaClY(layer.fMaxSigmaClY),
1636 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1637 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1642 //------------------------------------------------------------------------
1643 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1644 //--------------------------------------------------------------------
1645 // AliITSlayer destructor
1646 //--------------------------------------------------------------------
1647 delete [] fDetectors;
1648 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1649 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1650 fClusterWeight[i]=0;
1651 fClusterTracks[0][i]=-1;
1652 fClusterTracks[1][i]=-1;
1653 fClusterTracks[2][i]=-1;
1654 fClusterTracks[3][i]=-1;
1657 //------------------------------------------------------------------------
1658 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1659 //--------------------------------------------------------------------
1660 // This function removes loaded clusters
1661 //--------------------------------------------------------------------
1662 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1663 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1664 fClusterWeight[i]=0;
1665 fClusterTracks[0][i]=-1;
1666 fClusterTracks[1][i]=-1;
1667 fClusterTracks[2][i]=-1;
1668 fClusterTracks[3][i]=-1;
1674 //------------------------------------------------------------------------
1675 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1676 //--------------------------------------------------------------------
1677 // This function reset weights of the clusters
1678 //--------------------------------------------------------------------
1679 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1680 fClusterWeight[i]=0;
1681 fClusterTracks[0][i]=-1;
1682 fClusterTracks[1][i]=-1;
1683 fClusterTracks[2][i]=-1;
1684 fClusterTracks[3][i]=-1;
1686 for (Int_t i=0; i<fN;i++) {
1687 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1688 if (cl&&cl->IsUsed()) cl->Use();
1692 //------------------------------------------------------------------------
1693 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1694 //--------------------------------------------------------------------
1695 // This function calculates the road defined by the cluster density
1696 //--------------------------------------------------------------------
1698 for (Int_t i=0; i<fN; i++) {
1699 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1701 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1703 //------------------------------------------------------------------------
1704 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1705 //--------------------------------------------------------------------
1706 //This function adds a cluster to this layer
1707 //--------------------------------------------------------------------
1708 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1714 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1716 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1717 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1718 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1719 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1720 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1721 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1724 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1725 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1726 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1727 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1731 //------------------------------------------------------------------------
1732 void AliITStrackerMI::AliITSlayer::SortClusters()
1737 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1738 Float_t *z = new Float_t[fN];
1739 Int_t * index = new Int_t[fN];
1741 fMaxSigmaClY=0.; //AD
1742 fMaxSigmaClZ=0.; //AD
1744 for (Int_t i=0;i<fN;i++){
1745 z[i] = fClusters[i]->GetZ();
1746 // save largest errors in y and z for this layer
1747 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1748 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1750 TMath::Sort(fN,z,index,kFALSE);
1751 for (Int_t i=0;i<fN;i++){
1752 clusters[i] = fClusters[index[i]];
1755 for (Int_t i=0;i<fN;i++){
1756 fClusters[i] = clusters[i];
1757 fZ[i] = fClusters[i]->GetZ();
1758 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1759 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1760 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1770 for (Int_t i=0;i<fN;i++){
1771 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1772 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1773 fClusterIndex[i] = i;
1777 fDy5 = (fYB[1]-fYB[0])/5.;
1778 fDy10 = (fYB[1]-fYB[0])/10.;
1779 fDy20 = (fYB[1]-fYB[0])/20.;
1780 for (Int_t i=0;i<6;i++) fN5[i] =0;
1781 for (Int_t i=0;i<11;i++) fN10[i]=0;
1782 for (Int_t i=0;i<21;i++) fN20[i]=0;
1784 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;}
1785 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;}
1786 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;}
1789 for (Int_t i=0;i<fN;i++)
1790 for (Int_t irot=-1;irot<=1;irot++){
1791 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1793 for (Int_t slice=0; slice<6;slice++){
1794 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1795 fClusters5[slice][fN5[slice]] = fClusters[i];
1796 fY5[slice][fN5[slice]] = curY;
1797 fZ5[slice][fN5[slice]] = fZ[i];
1798 fClusterIndex5[slice][fN5[slice]]=i;
1803 for (Int_t slice=0; slice<11;slice++){
1804 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1805 fClusters10[slice][fN10[slice]] = fClusters[i];
1806 fY10[slice][fN10[slice]] = curY;
1807 fZ10[slice][fN10[slice]] = fZ[i];
1808 fClusterIndex10[slice][fN10[slice]]=i;
1813 for (Int_t slice=0; slice<21;slice++){
1814 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1815 fClusters20[slice][fN20[slice]] = fClusters[i];
1816 fY20[slice][fN20[slice]] = curY;
1817 fZ20[slice][fN20[slice]] = fZ[i];
1818 fClusterIndex20[slice][fN20[slice]]=i;
1825 // consistency check
1827 for (Int_t i=0;i<fN-1;i++){
1833 for (Int_t slice=0;slice<21;slice++)
1834 for (Int_t i=0;i<fN20[slice]-1;i++){
1835 if (fZ20[slice][i]>fZ20[slice][i+1]){
1842 //------------------------------------------------------------------------
1843 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1844 //--------------------------------------------------------------------
1845 // This function returns the index of the nearest cluster
1846 //--------------------------------------------------------------------
1849 if (fCurrentSlice<0) {
1858 if (ncl==0) return 0;
1859 Int_t b=0, e=ncl-1, m=(b+e)/2;
1860 for (; b<e; m=(b+e)/2) {
1861 // if (z > fClusters[m]->GetZ()) b=m+1;
1862 if (z > zcl[m]) b=m+1;
1867 //------------------------------------------------------------------------
1868 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 {
1869 //--------------------------------------------------------------------
1870 // This function computes the rectangular road for this track
1871 //--------------------------------------------------------------------
1874 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1875 // take into account the misalignment: propagate track to misaligned detector plane
1876 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1878 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1879 TMath::Sqrt(track->GetSigmaZ2() +
1880 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1881 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1882 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1883 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1884 TMath::Sqrt(track->GetSigmaY2() +
1885 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1886 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1887 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1889 // track at boundary between detectors, enlarge road
1890 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1891 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1892 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1893 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1894 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1895 Float_t tgl = TMath::Abs(track->GetTgl());
1896 if (tgl > 1.) tgl=1.;
1897 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1898 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1899 Float_t snp = TMath::Abs(track->GetSnp());
1900 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1901 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1904 // add to the road a term (up to 2-3 mm) to deal with misalignments
1905 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1906 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1908 Double_t r = fgLayers[ilayer].GetR();
1909 zmin = track->GetZ() - dz;
1910 zmax = track->GetZ() + dz;
1911 ymin = track->GetY() + r*det.GetPhi() - dy;
1912 ymax = track->GetY() + r*det.GetPhi() + dy;
1914 // bring track back to idead detector plane
1915 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1919 //------------------------------------------------------------------------
1920 void AliITStrackerMI::AliITSlayer::
1921 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1922 //--------------------------------------------------------------------
1923 // This function sets the "window"
1924 //--------------------------------------------------------------------
1926 Double_t circle=2*TMath::Pi()*fR;
1932 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1933 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1934 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1935 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1936 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1938 Float_t ymiddle = (fYmin+fYmax)*0.5;
1939 if (ymiddle<fYB[0]) {
1940 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1941 } else if (ymiddle>fYB[1]) {
1942 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1948 fClustersCs = fClusters;
1949 fClusterIndexCs = fClusterIndex;
1955 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1956 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1957 if (slice<0) slice=0;
1958 if (slice>20) slice=20;
1959 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1961 fCurrentSlice=slice;
1962 fClustersCs = fClusters20[fCurrentSlice];
1963 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1964 fYcs = fY20[fCurrentSlice];
1965 fZcs = fZ20[fCurrentSlice];
1966 fNcs = fN20[fCurrentSlice];
1971 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1972 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1973 if (slice<0) slice=0;
1974 if (slice>10) slice=10;
1975 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1977 fCurrentSlice=slice;
1978 fClustersCs = fClusters10[fCurrentSlice];
1979 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1980 fYcs = fY10[fCurrentSlice];
1981 fZcs = fZ10[fCurrentSlice];
1982 fNcs = fN10[fCurrentSlice];
1987 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1988 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1989 if (slice<0) slice=0;
1990 if (slice>5) slice=5;
1991 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1993 fCurrentSlice=slice;
1994 fClustersCs = fClusters5[fCurrentSlice];
1995 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1996 fYcs = fY5[fCurrentSlice];
1997 fZcs = fZ5[fCurrentSlice];
1998 fNcs = fN5[fCurrentSlice];
2002 fI = FindClusterIndex(fZmin);
2003 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2009 //------------------------------------------------------------------------
2010 Int_t AliITStrackerMI::AliITSlayer::
2011 FindDetectorIndex(Double_t phi, Double_t z) const {
2012 //--------------------------------------------------------------------
2013 //This function finds the detector crossed by the track
2014 //--------------------------------------------------------------------
2016 if (fZOffset<0) // old geometry
2017 dphi = -(phi-fPhiOffset);
2018 else // new geometry
2019 dphi = phi-fPhiOffset;
2022 if (dphi < 0) dphi += 2*TMath::Pi();
2023 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2024 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2025 if (np>=fNladders) np-=fNladders;
2026 if (np<0) np+=fNladders;
2029 Double_t dz=fZOffset-z;
2030 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2031 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2032 if (nz>=fNdetectors || nz<0) {
2033 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2037 // ad hoc correction for 3rd ladder of SDD inner layer,
2038 // which is reversed (rotated by pi around local y)
2039 // this correction is OK only from AliITSv11Hybrid onwards
2040 if (GetR()>12. && GetR()<20.) { // SDD inner
2041 if(np==2) { // 3rd ladder
2042 Double_t posMod252[3];
2043 AliITSgeomTGeo::GetTranslation(252,posMod252);
2044 // check the Z coordinate of Mod 252: if negative
2045 // (old SDD geometry in AliITSv11Hybrid)
2046 // the swap of numeration whould be applied
2047 if(posMod252[2]<0.){
2048 nz = (fNdetectors-1) - nz;
2052 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
2055 return np*fNdetectors + nz;
2057 //------------------------------------------------------------------------
2058 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2060 //--------------------------------------------------------------------
2061 // This function returns clusters within the "window"
2062 //--------------------------------------------------------------------
2064 if (fCurrentSlice<0) {
2065 Double_t rpi2 = 2.*fR*TMath::Pi();
2066 for (Int_t i=fI; i<fImax; i++) {
2069 if (fYmax<y) y -= rpi2;
2070 if (fYmin>y) y += rpi2;
2071 if (y<fYmin) continue;
2072 if (y>fYmax) continue;
2074 // skip clusters that are in "extended" road but they
2075 // 3sigma error does not touch the original road
2076 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2077 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2079 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2082 return fClusters[i];
2085 for (Int_t i=fI; i<fImax; i++) {
2086 if (fYcs[i]<fYmin) continue;
2087 if (fYcs[i]>fYmax) continue;
2088 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2089 ci=fClusterIndexCs[i];
2091 return fClustersCs[i];
2096 //------------------------------------------------------------------------
2097 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2099 //--------------------------------------------------------------------
2100 // This function returns the layer thickness at this point (units X0)
2101 //--------------------------------------------------------------------
2103 x0=AliITSRecoParam::GetX0Air();
2104 if (43<fR&&fR<45) { //SSD2
2107 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2108 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2109 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2110 for (Int_t i=0; i<12; i++) {
2111 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2112 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2116 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2117 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2121 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2122 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2125 if (37<fR&&fR<41) { //SSD1
2128 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2129 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2130 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2131 for (Int_t i=0; i<11; i++) {
2132 if (TMath::Abs(z-3.9*i)<0.15) {
2133 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2137 if (TMath::Abs(z+3.9*i)<0.15) {
2138 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2142 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2143 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2146 if (13<fR&&fR<26) { //SDD
2149 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2151 if (TMath::Abs(y-1.80)<0.55) {
2153 for (Int_t j=0; j<20; j++) {
2154 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2155 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2158 if (TMath::Abs(y+1.80)<0.55) {
2160 for (Int_t j=0; j<20; j++) {
2161 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2162 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2166 for (Int_t i=0; i<4; i++) {
2167 if (TMath::Abs(z-7.3*i)<0.60) {
2169 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2172 if (TMath::Abs(z+7.3*i)<0.60) {
2174 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2179 if (6<fR&&fR<8) { //SPD2
2180 Double_t dd=0.0063; x0=21.5;
2182 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2183 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2185 if (3<fR&&fR<5) { //SPD1
2186 Double_t dd=0.0063; x0=21.5;
2188 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2189 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2194 //------------------------------------------------------------------------
2195 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2197 fRmisal(det.fRmisal),
2199 fSinPhi(det.fSinPhi),
2200 fCosPhi(det.fCosPhi),
2206 fNChips(det.fNChips),
2207 fChipIsBad(det.fChipIsBad)
2211 //------------------------------------------------------------------------
2212 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2213 const AliITSDetTypeRec *detTypeRec)
2215 //--------------------------------------------------------------------
2216 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2217 //--------------------------------------------------------------------
2219 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2220 // while in the tracker they start from 0 for each layer
2221 for(Int_t il=0; il<ilayer; il++)
2222 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2225 if (ilayer==0 || ilayer==1) { // ---------- SPD
2227 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2229 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2232 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2236 // Get calibration from AliITSDetTypeRec
2237 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2238 calib->SetModuleIndex(idet);
2239 AliITSCalibration *calibSPDdead = 0;
2240 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2241 if (calib->IsBad() ||
2242 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2245 // printf("lay %d bad %d\n",ilayer,idet);
2248 // Get segmentation from AliITSDetTypeRec
2249 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2251 // Read info about bad chips
2252 fNChips = segm->GetMaximumChipIndex()+1;
2253 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2254 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2255 fChipIsBad = new Bool_t[fNChips];
2256 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2257 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2258 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2259 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2264 //------------------------------------------------------------------------
2265 Double_t AliITStrackerMI::GetEffectiveThickness()
2267 //--------------------------------------------------------------------
2268 // Returns the thickness between the current layer and the vertex (units X0)
2269 //--------------------------------------------------------------------
2272 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2273 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2274 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2278 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2279 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2283 Double_t xn=fgLayers[fI].GetR();
2284 for (Int_t i=0; i<fI; i++) {
2285 Double_t xi=fgLayers[i].GetR();
2286 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2292 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2293 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2296 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2297 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2301 //------------------------------------------------------------------------
2302 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2303 //-------------------------------------------------------------------
2304 // This function returns number of clusters within the "window"
2305 //--------------------------------------------------------------------
2307 for (Int_t i=fI; i<fN; i++) {
2308 const AliITSRecPoint *c=fClusters[i];
2309 if (c->GetZ() > fZmax) break;
2310 if (c->IsUsed()) continue;
2311 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2312 Double_t y=fR*det.GetPhi() + c->GetY();
2314 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2315 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2317 if (y<fYmin) continue;
2318 if (y>fYmax) continue;
2323 //------------------------------------------------------------------------
2324 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2325 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2327 //--------------------------------------------------------------------
2328 // This function refits the track "track" at the position "x" using
2329 // the clusters from "clusters"
2330 // If "extra"==kTRUE,
2331 // the clusters from overlapped modules get attached to "track"
2332 // If "planeff"==kTRUE,
2333 // special approach for plane efficiency evaluation is applyed
2334 //--------------------------------------------------------------------
2336 Int_t index[AliITSgeomTGeo::kNLayers];
2338 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2339 Int_t nc=clusters->GetNumberOfClusters();
2340 for (k=0; k<nc; k++) {
2341 Int_t idx=clusters->GetClusterIndex(k);
2342 Int_t ilayer=(idx&0xf0000000)>>28;
2346 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2348 //------------------------------------------------------------------------
2349 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2350 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2352 //--------------------------------------------------------------------
2353 // This function refits the track "track" at the position "x" using
2354 // the clusters from array
2355 // If "extra"==kTRUE,
2356 // the clusters from overlapped modules get attached to "track"
2357 // If "planeff"==kTRUE,
2358 // special approach for plane efficiency evaluation is applyed
2359 //--------------------------------------------------------------------
2360 Int_t index[AliITSgeomTGeo::kNLayers];
2362 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2364 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2365 index[k]=clusters[k];
2368 // special for cosmics and TPC prolonged tracks:
2369 // propagate to the innermost of:
2370 // - innermost layer crossed by the track
2371 // - innermost layer where a cluster was associated to the track
2372 static AliITSRecoParam *repa = NULL;
2374 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2376 repa = AliITSRecoParam::GetHighFluxParam();
2377 AliWarning("Using default AliITSRecoParam class");
2380 Int_t evsp=repa->GetEventSpecie();
2382 if(track->GetESDtrack()) trStatus=track->GetStatus();
2383 Int_t innermostlayer=0;
2384 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2386 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2387 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2388 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2389 index[innermostlayer] >= 0 ) break;
2392 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2395 Int_t modstatus=1; // found
2397 Int_t from, to, step;
2398 if (xx > track->GetX()) {
2399 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2402 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2405 TString dir = (step>0 ? "outward" : "inward");
2407 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2408 AliITSlayer &layer=fgLayers[ilayer];
2409 Double_t r=layer.GetR();
2411 if (step<0 && xx>r) break;
2413 // material between SSD and SDD, SDD and SPD
2414 Double_t hI=ilayer-0.5*step;
2415 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2416 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2417 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2418 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2421 Double_t oldGlobXYZ[3];
2422 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2424 // continue if we are already beyond this layer
2425 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2426 if(step>0 && oldGlobR > r) continue; // going outward
2427 if(step<0 && oldGlobR < r) continue; // going inward
2430 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2432 Int_t idet=layer.FindDetectorIndex(phi,z);
2434 // check if we allow a prolongation without point for large-eta tracks
2435 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2437 modstatus = 4; // out in z
2438 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2439 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2442 // apply correction for material of the current layer
2443 // add time if going outward
2444 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2448 if (idet<0) return kFALSE;
2451 const AliITSdetector &det=layer.GetDetector(idet);
2452 // only for ITS-SA tracks refit
2453 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2455 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2457 track->SetDetectorIndex(idet);
2458 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2460 Double_t dz,zmin,zmax,dy,ymin,ymax;
2462 const AliITSRecPoint *clAcc=0;
2463 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2465 Int_t idx=index[ilayer];
2466 if (idx>=0) { // cluster in this layer
2467 modstatus = 6; // no refit
2468 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2470 if (idet != cl->GetDetectorIndex()) {
2471 idet=cl->GetDetectorIndex();
2472 const AliITSdetector &detc=layer.GetDetector(idet);
2473 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2474 track->SetDetectorIndex(idet);
2475 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2477 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2478 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2482 modstatus = 1; // found
2487 } else { // no cluster in this layer
2489 modstatus = 3; // skipped
2490 // Plane Eff determination:
2491 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2492 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2493 UseTrackForPlaneEff(track,ilayer);
2496 modstatus = 5; // no cls in road
2498 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2499 dz = 0.5*(zmax-zmin);
2500 dy = 0.5*(ymax-ymin);
2501 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2502 if (dead==1) modstatus = 7; // holes in z in SPD
2503 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2508 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2509 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2511 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2514 if (extra && clAcc) { // search for extra clusters in overlapped modules
2515 AliITStrackV2 tmp(*track);
2516 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2517 layer.SelectClusters(zmin,zmax,ymin,ymax);
2519 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2521 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2522 Double_t tolerance=0.1;
2523 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2524 // only clusters in another module! (overlaps)
2525 idetExtra = clExtra->GetDetectorIndex();
2526 if (idet == idetExtra) continue;
2528 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2530 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2531 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2532 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2533 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2535 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2536 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2539 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2540 track->SetExtraModule(ilayer,idetExtra);
2542 } // end search for extra clusters in overlapped modules
2544 // Correct for material of the current layer
2546 // add time if going outward
2547 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2548 track->SetCheckInvariant(kTRUE);
2549 } // end loop on layers
2551 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2555 //------------------------------------------------------------------------
2556 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2559 // calculate normalized chi2
2560 // return NormalizedChi2(track,0);
2563 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2564 // track->fdEdxMismatch=0;
2565 Float_t dedxmismatch =0;
2566 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2568 for (Int_t i = 0;i<6;i++){
2569 if (track->GetClIndex(i)>=0){
2570 Float_t cerry, cerrz;
2571 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2573 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2576 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2577 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2578 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2580 cchi2+=(0.5-ratio)*10.;
2581 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2582 dedxmismatch+=(0.5-ratio)*10.;
2586 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2587 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2588 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2589 if (i<2) chi2+=2*cl->GetDeltaProbability();
2595 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2596 track->SetdEdxMismatch(dedxmismatch);
2600 for (Int_t i = 0;i<4;i++){
2601 if (track->GetClIndex(i)>=0){
2602 Float_t cerry, cerrz;
2603 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2604 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2607 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2608 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2612 for (Int_t i = 4;i<6;i++){
2613 if (track->GetClIndex(i)>=0){
2614 Float_t cerry, cerrz;
2615 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2616 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2619 Float_t cerryb, cerrzb;
2620 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2621 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2624 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2625 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2630 if (track->GetESDtrack()->GetTPCsignal()>85){
2631 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2633 chi2+=(0.5-ratio)*5.;
2636 chi2+=(ratio-2.0)*3;
2640 Double_t match = TMath::Sqrt(track->GetChi22());
2641 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2642 if (!track->GetConstrain()) {
2643 if (track->GetNumberOfClusters()>2) {
2644 match/=track->GetNumberOfClusters()-2.;
2649 if (match<0) match=0;
2651 // penalty factor for missing points (NDeadZone>0), but no penalty
2652 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2653 // or there is a dead from OCDB)
2654 Float_t deadzonefactor = 0.;
2655 if (track->GetNDeadZone()>0.) {
2656 Int_t sumDeadZoneProbability=0;
2657 for(Int_t ilay=0;ilay<6;ilay++) {
2658 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2660 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2661 if(nDeadZoneWithProbNot1>0) {
2662 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2663 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2664 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2666 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2667 deadzonefactor = 3.*(1.1-deadZoneProbability);
2671 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2672 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2673 1./(1.+track->GetNSkipped()));
2674 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2675 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2678 //------------------------------------------------------------------------
2679 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2682 // return matching chi2 between two tracks
2683 Double_t largeChi2=1000.;
2685 AliITStrackMI track3(*track2);
2686 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2688 vec(0,0)=track1->GetY() - track3.GetY();
2689 vec(1,0)=track1->GetZ() - track3.GetZ();
2690 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2691 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2692 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2695 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2696 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2697 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2698 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2699 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2701 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2702 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2703 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2704 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2706 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2707 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2708 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2710 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2711 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2713 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2716 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2717 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2720 //------------------------------------------------------------------------
2721 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2724 // return probability that given point (characterized by z position and error)
2725 // is in SPD dead zone
2726 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2728 Double_t probability = 0.;
2729 Double_t nearestz = 0.,distz=0.;
2730 Int_t nearestzone = -1;
2731 Double_t mindistz = 1000.;
2733 // find closest dead zone
2734 for (Int_t i=0; i<3; i++) {
2735 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2736 if (distz<mindistz) {
2738 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2743 // too far from dead zone
2744 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2747 Double_t zmin, zmax;
2748 if (nearestzone==0) { // dead zone at z = -7
2749 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2750 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2751 } else if (nearestzone==1) { // dead zone at z = 0
2752 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2753 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2754 } else if (nearestzone==2) { // dead zone at z = +7
2755 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2756 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2761 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2763 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2764 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2765 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2768 //------------------------------------------------------------------------
2769 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2772 // calculate normalized chi2
2774 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2776 for (Int_t i = 0;i<6;i++){
2777 if (TMath::Abs(track->GetDy(i))>0){
2778 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2779 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2782 else{chi2[i]=10000;}
2785 TMath::Sort(6,chi2,index,kFALSE);
2786 Float_t max = float(ncl)*fac-1.;
2787 Float_t sumchi=0, sumweight=0;
2788 for (Int_t i=0;i<max+1;i++){
2789 Float_t weight = (i<max)?1.:(max+1.-i);
2790 sumchi+=weight*chi2[index[i]];
2793 Double_t normchi2 = sumchi/sumweight;
2796 //------------------------------------------------------------------------
2797 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2800 // calculate normalized chi2
2801 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2804 for (Int_t i=0;i<6;i++){
2805 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2806 Double_t sy1 = forwardtrack->GetSigmaY(i);
2807 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2808 Double_t sy2 = backtrack->GetSigmaY(i);
2809 Double_t sz2 = backtrack->GetSigmaZ(i);
2810 if (i<2){ sy2=1000.;sz2=1000;}
2812 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2813 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2815 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2816 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2818 res+= nz0*nz0+ny0*ny0;
2821 if (npoints>1) return
2822 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2823 //2*forwardtrack->fNUsed+
2824 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2825 1./(1.+forwardtrack->GetNSkipped()));
2828 //------------------------------------------------------------------------
2829 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2830 //--------------------------------------------------------------------
2831 // Return pointer to a given cluster
2832 //--------------------------------------------------------------------
2833 Int_t l=(index & 0xf0000000) >> 28;
2834 Int_t c=(index & 0x0fffffff) >> 00;
2835 return fgLayers[l].GetWeight(c);
2837 //------------------------------------------------------------------------
2838 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2840 //---------------------------------------------
2841 // register track to the list
2843 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2846 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2847 Int_t index = track->GetClusterIndex(icluster);
2848 Int_t l=(index & 0xf0000000) >> 28;
2849 Int_t c=(index & 0x0fffffff) >> 00;
2850 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2851 for (Int_t itrack=0;itrack<4;itrack++){
2852 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2853 fgLayers[l].SetClusterTracks(itrack,c,id);
2859 //------------------------------------------------------------------------
2860 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2862 //---------------------------------------------
2863 // unregister track from the list
2864 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2865 Int_t index = track->GetClusterIndex(icluster);
2866 Int_t l=(index & 0xf0000000) >> 28;
2867 Int_t c=(index & 0x0fffffff) >> 00;
2868 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2869 for (Int_t itrack=0;itrack<4;itrack++){
2870 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2871 fgLayers[l].SetClusterTracks(itrack,c,-1);
2876 //------------------------------------------------------------------------
2877 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2879 //-------------------------------------------------------------
2880 //get number of shared clusters
2881 //-------------------------------------------------------------
2883 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2884 // mean number of clusters
2885 Float_t *ny = GetNy(id), *nz = GetNz(id);
2888 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2889 Int_t index = track->GetClusterIndex(icluster);
2890 Int_t l=(index & 0xf0000000) >> 28;
2891 Int_t c=(index & 0x0fffffff) >> 00;
2892 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2893 // if (ny[l]<1.e-13){
2894 // printf("problem\n");
2896 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2900 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2901 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2902 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2904 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2907 deltan = (cl->GetNz()-nz[l]);
2909 if (deltan>2.0) continue; // extended - highly probable shared cluster
2910 weight = 2./TMath::Max(3.+deltan,2.);
2912 for (Int_t itrack=0;itrack<4;itrack++){
2913 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2915 clist[l] = (AliITSRecPoint*)GetCluster(index);
2916 track->SetSharedWeight(l,weight);
2922 track->SetNUsed(shared);
2925 //------------------------------------------------------------------------
2926 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2929 // find first shared track
2931 // mean number of clusters
2932 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2934 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2935 Int_t sharedtrack=100000;
2936 Int_t tracks[24],trackindex=0;
2937 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2939 for (Int_t icluster=0;icluster<6;icluster++){
2940 if (clusterlist[icluster]<0) continue;
2941 Int_t index = clusterlist[icluster];
2942 Int_t l=(index & 0xf0000000) >> 28;
2943 Int_t c=(index & 0x0fffffff) >> 00;
2944 // if (ny[l]<1.e-13){
2945 // printf("problem\n");
2947 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2948 //if (l>3) continue;
2949 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2952 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2953 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2954 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2956 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2959 deltan = (cl->GetNz()-nz[l]);
2961 if (deltan>2.0) continue; // extended - highly probable shared cluster
2963 for (Int_t itrack=3;itrack>=0;itrack--){
2964 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2965 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2966 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2971 if (trackindex==0) return -1;
2973 sharedtrack = tracks[0];
2975 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2978 Int_t tracks2[24], cluster[24];
2979 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2982 for (Int_t i=0;i<trackindex;i++){
2983 if (tracks[i]<0) continue;
2984 tracks2[index] = tracks[i];
2986 for (Int_t j=i+1;j<trackindex;j++){
2987 if (tracks[j]<0) continue;
2988 if (tracks[j]==tracks[i]){
2996 for (Int_t i=0;i<index;i++){
2997 if (cluster[index]>max) {
2998 sharedtrack=tracks2[index];
3005 if (sharedtrack>=100000) return -1;
3007 // make list of overlaps
3009 for (Int_t icluster=0;icluster<6;icluster++){
3010 if (clusterlist[icluster]<0) continue;
3011 Int_t index = clusterlist[icluster];
3012 Int_t l=(index & 0xf0000000) >> 28;
3013 Int_t c=(index & 0x0fffffff) >> 00;
3014 if (c>fgLayers[l].GetNumberOfClusters()) continue;
3015 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3017 if (cl->GetNy()>2) continue;
3018 if (cl->GetNz()>2) continue;
3021 if (cl->GetNy()>3) continue;
3022 if (cl->GetNz()>3) continue;
3025 for (Int_t itrack=3;itrack>=0;itrack--){
3026 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3027 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3035 //------------------------------------------------------------------------
3036 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
3038 // try to find track hypothesys without conflicts
3039 // with minimal chi2;
3040 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3041 Int_t entries1 = arr1->GetEntriesFast();
3042 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3043 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3044 Int_t entries2 = arr2->GetEntriesFast();
3045 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3047 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3048 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3049 if (track10->Pt()>0.5+track20->Pt()) return track10;
3051 for (Int_t itrack=0;itrack<entries1;itrack++){
3052 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3053 UnRegisterClusterTracks(track,trackID1);
3056 for (Int_t itrack=0;itrack<entries2;itrack++){
3057 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3058 UnRegisterClusterTracks(track,trackID2);
3062 Float_t maxconflicts=6;
3063 Double_t maxchi2 =1000.;
3065 // get weight of hypothesys - using best hypothesy
3068 Int_t list1[6],list2[6];
3069 AliITSRecPoint *clist1[6], *clist2[6] ;
3070 RegisterClusterTracks(track10,trackID1);
3071 RegisterClusterTracks(track20,trackID2);
3072 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3073 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3074 UnRegisterClusterTracks(track10,trackID1);
3075 UnRegisterClusterTracks(track20,trackID2);
3078 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3079 Float_t nerry[6],nerrz[6];
3080 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3081 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3082 for (Int_t i=0;i<6;i++){
3083 if ( (erry1[i]>0) && (erry2[i]>0)) {
3084 nerry[i] = TMath::Min(erry1[i],erry2[i]);
3085 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3087 nerry[i] = TMath::Max(erry1[i],erry2[i]);
3088 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3090 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3091 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3092 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3095 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3096 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3097 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3105 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3106 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3107 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3108 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3110 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3111 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3112 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3114 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3115 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3116 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3119 Double_t sumw = w1+w2;
3123 w1 = (d2+0.5)/(d1+d2+1.);
3124 w2 = (d1+0.5)/(d1+d2+1.);
3126 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3127 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3129 // get pair of "best" hypothesys
3131 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3132 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3134 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3135 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3136 //if (track1->fFakeRatio>0) continue;
3137 RegisterClusterTracks(track1,trackID1);
3138 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3139 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3141 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3142 //if (track2->fFakeRatio>0) continue;
3144 RegisterClusterTracks(track2,trackID2);
3145 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3146 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3147 UnRegisterClusterTracks(track2,trackID2);
3149 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3150 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3151 if (nskipped>0.5) continue;
3153 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3154 if (conflict1+1<cconflict1) continue;
3155 if (conflict2+1<cconflict2) continue;
3159 for (Int_t i=0;i<6;i++){
3161 Float_t c1 =0.; // conflict coeficients
3163 if (clist1[i]&&clist2[i]){
3166 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3169 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3171 c1 = 2./TMath::Max(3.+deltan,2.);
3172 c2 = 2./TMath::Max(3.+deltan,2.);
3178 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3181 deltan = (clist1[i]->GetNz()-nz1[i]);
3183 c1 = 2./TMath::Max(3.+deltan,2.);
3190 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3193 deltan = (clist2[i]->GetNz()-nz2[i]);
3195 c2 = 2./TMath::Max(3.+deltan,2.);
3201 if (TMath::Abs(track1->GetDy(i))>0.) {
3202 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3203 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3204 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3205 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3207 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3210 if (TMath::Abs(track2->GetDy(i))>0.) {
3211 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3212 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3213 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3214 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3217 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3219 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3220 if (chi21>0) sum+=w1;
3221 if (chi22>0) sum+=w2;
3224 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3225 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3226 Double_t normchi2 = 2*conflict+sumchi2/sum;
3227 if ( normchi2 <maxchi2 ){
3230 maxconflicts = conflict;
3234 UnRegisterClusterTracks(track1,trackID1);
3237 // if (maxconflicts<4 && maxchi2<th0){
3238 if (maxchi2<th0*2.){
3239 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3240 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3241 track1->SetChi2MIP(5,maxconflicts);
3242 track1->SetChi2MIP(6,maxchi2);
3243 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3244 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3245 track1->SetChi2MIP(8,index1);
3246 fBestTrackIndex[trackID1] =index1;
3247 UpdateESDtrack(track1, AliESDtrack::kITSin);
3249 else if (track10->GetChi2MIP(0)<th1){
3250 track10->SetChi2MIP(5,maxconflicts);
3251 track10->SetChi2MIP(6,maxchi2);
3252 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3253 UpdateESDtrack(track10,AliESDtrack::kITSin);
3256 for (Int_t itrack=0;itrack<entries1;itrack++){
3257 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3258 UnRegisterClusterTracks(track,trackID1);
3261 for (Int_t itrack=0;itrack<entries2;itrack++){
3262 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3263 UnRegisterClusterTracks(track,trackID2);
3266 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3267 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3268 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3269 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3270 RegisterClusterTracks(track10,trackID1);
3272 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3273 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3274 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3275 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3276 RegisterClusterTracks(track20,trackID2);
3281 //------------------------------------------------------------------------
3282 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3283 //--------------------------------------------------------------------
3284 // This function marks clusters assigned to the track
3285 //--------------------------------------------------------------------
3286 AliTracker::UseClusters(t,from);
3288 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3289 //if (c->GetQ()>2) c->Use();
3290 if (c->GetSigmaZ2()>0.1) c->Use();
3291 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3292 //if (c->GetQ()>2) c->Use();
3293 if (c->GetSigmaZ2()>0.1) c->Use();
3296 //------------------------------------------------------------------------
3297 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3299 //------------------------------------------------------------------
3300 // add track to the list of hypothesys
3301 //------------------------------------------------------------------
3303 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3304 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3306 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3308 array = new TObjArray(10);
3309 fTrackHypothesys.AddAt(array,esdindex);
3311 array->AddLast(track);
3313 //------------------------------------------------------------------------
3314 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3316 //-------------------------------------------------------------------
3317 // compress array of track hypothesys
3318 // keep only maxsize best hypothesys
3319 //-------------------------------------------------------------------
3320 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3321 if (! (fTrackHypothesys.At(esdindex)) ) return;
3322 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3323 Int_t entries = array->GetEntriesFast();
3325 //- find preliminary besttrack as a reference
3326 Float_t minchi2=10000;
3328 AliITStrackMI * besttrack=0;
3329 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3330 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3331 if (!track) continue;
3332 Float_t chi2 = NormalizedChi2(track,0);
3334 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3335 track->SetLabel(tpcLabel);
3337 track->SetFakeRatio(1.);
3338 CookLabel(track,0.); //For comparison only
3340 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3341 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3342 if (track->GetNumberOfClusters()<maxn) continue;
3343 maxn = track->GetNumberOfClusters();
3350 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3351 delete array->RemoveAt(itrack);
3355 if (!besttrack) return;
3358 //take errors of best track as a reference
3359 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3360 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3361 for (Int_t j=0;j<6;j++) {
3362 if (besttrack->GetClIndex(j)>=0){
3363 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3364 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3365 ny[j] = besttrack->GetNy(j);
3366 nz[j] = besttrack->GetNz(j);
3370 // calculate normalized chi2
3372 Float_t * chi2 = new Float_t[entries];
3373 Int_t * index = new Int_t[entries];
3374 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3375 for (Int_t itrack=0;itrack<entries;itrack++){
3376 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3378 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3379 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3380 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3381 chi2[itrack] = track->GetChi2MIP(0);
3383 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3384 delete array->RemoveAt(itrack);
3390 TMath::Sort(entries,chi2,index,kFALSE);
3391 besttrack = (AliITStrackMI*)array->At(index[0]);
3392 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3393 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3394 for (Int_t j=0;j<6;j++){
3395 if (besttrack->GetClIndex(j)>=0){
3396 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3397 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3398 ny[j] = besttrack->GetNy(j);
3399 nz[j] = besttrack->GetNz(j);
3404 // calculate one more time with updated normalized errors
3405 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3406 for (Int_t itrack=0;itrack<entries;itrack++){
3407 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3409 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3410 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3411 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3412 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3415 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3416 delete array->RemoveAt(itrack);
3421 entries = array->GetEntriesFast();
3425 TObjArray * newarray = new TObjArray();
3426 TMath::Sort(entries,chi2,index,kFALSE);
3427 besttrack = (AliITStrackMI*)array->At(index[0]);
3429 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3431 for (Int_t j=0;j<6;j++){
3432 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3433 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3434 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3435 ny[j] = besttrack->GetNy(j);
3436 nz[j] = besttrack->GetNz(j);
3439 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3440 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3441 Float_t minn = besttrack->GetNumberOfClusters()-3;
3443 for (Int_t i=0;i<entries;i++){
3444 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3445 if (!track) continue;
3446 if (accepted>maxcut) break;
3447 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3448 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3449 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3450 delete array->RemoveAt(index[i]);
3454 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3455 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3456 if (!shortbest) accepted++;
3458 newarray->AddLast(array->RemoveAt(index[i]));
3459 for (Int_t j=0;j<6;j++){
3461 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3462 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3463 ny[j] = track->GetNy(j);
3464 nz[j] = track->GetNz(j);
3469 delete array->RemoveAt(index[i]);
3473 delete fTrackHypothesys.RemoveAt(esdindex);
3474 fTrackHypothesys.AddAt(newarray,esdindex);
3478 delete fTrackHypothesys.RemoveAt(esdindex);
3484 //------------------------------------------------------------------------
3485 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3487 //-------------------------------------------------------------
3488 // try to find best hypothesy
3489 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3490 //-------------------------------------------------------------
3491 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3492 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3493 if (!array) return 0;
3494 Int_t entries = array->GetEntriesFast();
3495 if (!entries) return 0;
3496 Float_t minchi2 = 100000;
3497 AliITStrackMI * besttrack=0;
3499 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3500 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3501 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3502 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3504 for (Int_t i=0;i<entries;i++){
3505 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3506 if (!track) continue;
3507 Float_t sigmarfi,sigmaz;
3508 GetDCASigma(track,sigmarfi,sigmaz);
3509 track->SetDnorm(0,sigmarfi);
3510 track->SetDnorm(1,sigmaz);
3512 track->SetChi2MIP(1,1000000);
3513 track->SetChi2MIP(2,1000000);
3514 track->SetChi2MIP(3,1000000);
3517 backtrack = new(backtrack) AliITStrackMI(*track);
3518 if (track->GetConstrain()) {
3519 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3520 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3521 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3523 backtrack->ResetCovariance(10.);
3525 backtrack->ResetCovariance(10.);
3527 backtrack->ResetClusters();
3529 Double_t x = original->GetX();
3530 if (!RefitAt(x,backtrack,track)) continue;
3532 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3533 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3534 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3535 track->SetChi22(GetMatchingChi2(backtrack,original));
3537 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3538 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3539 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3542 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3544 //forward track - without constraint
3545 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3546 forwardtrack->ResetClusters();
3548 RefitAt(x,forwardtrack,track);
3549 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3550 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3551 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3553 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3554 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3555 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3556 forwardtrack->SetD(0,track->GetD(0));
3557 forwardtrack->SetD(1,track->GetD(1));
3560 AliITSRecPoint* clist[6];
3561 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3562 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3565 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3566 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3567 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3568 track->SetChi2MIP(3,1000);
3571 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3573 for (Int_t ichi=0;ichi<5;ichi++){
3574 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3576 if (chi2 < minchi2){
3577 //besttrack = new AliITStrackMI(*forwardtrack);
3579 besttrack->SetLabel(track->GetLabel());
3580 besttrack->SetFakeRatio(track->GetFakeRatio());
3582 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3583 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3584 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3588 delete forwardtrack;
3590 if (!besttrack) return 0;
3593 for (Int_t i=0;i<entries;i++){
3594 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3596 if (!track) continue;
3598 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3599 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3600 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3601 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3602 delete array->RemoveAt(i);
3612 SortTrackHypothesys(esdindex,checkmax,1);
3614 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3615 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3616 besttrack = (AliITStrackMI*)array->At(0);
3617 if (!besttrack) return 0;
3618 besttrack->SetChi2MIP(8,0);
3619 fBestTrackIndex[esdindex]=0;
3620 entries = array->GetEntriesFast();
3621 AliITStrackMI *longtrack =0;
3623 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3624 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3625 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3626 if (!track->GetConstrain()) continue;
3627 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3628 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3629 if (track->GetChi2MIP(0)>4.) continue;
3630 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3633 //if (longtrack) besttrack=longtrack;
3636 AliITSRecPoint * clist[6];
3637 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3638 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3639 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3640 RegisterClusterTracks(besttrack,esdindex);
3647 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3648 if (sharedtrack>=0){
3650 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3652 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3658 if (shared>2.5) return 0;
3659 if (shared>1.0) return besttrack;
3661 // Don't sign clusters if not gold track
3663 if (!besttrack->IsGoldPrimary()) return besttrack;
3664 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3666 if (fConstraint[fPass]){
3670 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3671 for (Int_t i=0;i<6;i++){
3672 Int_t index = besttrack->GetClIndex(i);
3673 if (index<0) continue;
3674 Int_t ilayer = (index & 0xf0000000) >> 28;
3675 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3676 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3678 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3679 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3680 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3681 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3682 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3683 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3685 Bool_t cansign = kTRUE;
3686 for (Int_t itrack=0;itrack<entries; itrack++){
3687 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3688 if (!track) continue;
3689 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3690 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3696 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3697 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3698 if (!c->IsUsed()) c->Use();
3704 //------------------------------------------------------------------------
3705 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3708 // get "best" hypothesys
3711 Int_t nentries = itsTracks.GetEntriesFast();
3712 for (Int_t i=0;i<nentries;i++){
3713 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3714 if (!track) continue;
3715 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3716 if (!array) continue;
3717 if (array->GetEntriesFast()<=0) continue;
3719 AliITStrackMI* longtrack=0;
3721 Float_t maxchi2=1000;
3722 for (Int_t j=0;j<array->GetEntriesFast();j++){
3723 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3724 if (!trackHyp) continue;
3725 if (trackHyp->GetGoldV0()) {
3726 longtrack = trackHyp; //gold V0 track taken
3729 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3730 Float_t chi2 = trackHyp->GetChi2MIP(0);
3732 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3734 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3736 if (chi2 > maxchi2) continue;
3737 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3744 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3745 if (!longtrack) {longtrack = besttrack;}
3746 else besttrack= longtrack;
3750 AliITSRecPoint * clist[6];
3751 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3753 track->SetNUsed(shared);
3754 track->SetNSkipped(besttrack->GetNSkipped());
3755 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3757 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3761 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3762 //if (sharedtrack==-1) sharedtrack=0;
3763 if (sharedtrack>=0) {
3764 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3767 if (besttrack&&fAfterV0) {
3768 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3771 if (fConstraint[fPass]) UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3772 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3773 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3774 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3780 //------------------------------------------------------------------------
3781 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3782 //--------------------------------------------------------------------
3783 //This function "cooks" a track label. If label<0, this track is fake.
3784 //--------------------------------------------------------------------
3787 if (track->GetESDtrack()){
3788 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3789 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3790 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3792 track->SetChi2MIP(9,0);
3794 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3795 Int_t cindex = track->GetClusterIndex(i);
3796 Int_t l=(cindex & 0xf0000000) >> 28;
3797 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3799 for (Int_t ind=0;ind<3;ind++){
3800 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3801 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3803 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3806 Int_t nclusters = track->GetNumberOfClusters();
3807 if (nclusters > 0) //PH Some tracks don't have any cluster
3808 track->SetFakeRatio(double(nwrong)/double(nclusters));
3809 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3810 track->SetLabel(-tpcLabel);
3812 track->SetLabel(tpcLabel);
3814 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3817 //------------------------------------------------------------------------
3818 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
3820 // Fill the dE/dx in this track
3822 track->SetChi2MIP(9,0);
3823 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3824 Int_t cindex = track->GetClusterIndex(i);
3825 Int_t l=(cindex & 0xf0000000) >> 28;
3826 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3827 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3829 for (Int_t ind=0;ind<3;ind++){
3830 if (cl->GetLabel(ind)==lab) isWrong=0;
3832 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3836 track->CookdEdx(low,up);
3838 //------------------------------------------------------------------------
3839 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3841 // Create some arrays
3843 if (fCoefficients) delete []fCoefficients;
3844 fCoefficients = new Float_t[ntracks*48];
3845 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3847 //------------------------------------------------------------------------
3848 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3851 // Compute predicted chi2
3853 // Take into account the mis-alignment (bring track to cluster plane)
3854 Double_t xTrOrig=track->GetX();
3855 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3856 Float_t erry,errz,covyz;
3857 Float_t theta = track->GetTgl();
3858 Float_t phi = track->GetSnp();
3859 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3860 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3861 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()));
3862 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()));
3863 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3864 // Bring the track back to detector plane in ideal geometry
3865 // [mis-alignment will be accounted for in UpdateMI()]
3866 if (!track->Propagate(xTrOrig)) return 1000.;
3868 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3869 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3871 chi2+=0.5*TMath::Min(delta/2,2.);
3872 chi2+=2.*cluster->GetDeltaProbability();
3875 track->SetNy(layer,ny);
3876 track->SetNz(layer,nz);
3877 track->SetSigmaY(layer,erry);
3878 track->SetSigmaZ(layer, errz);
3879 track->SetSigmaYZ(layer,covyz);
3880 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3881 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3885 //------------------------------------------------------------------------
3886 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3891 Int_t layer = (index & 0xf0000000) >> 28;
3892 track->SetClIndex(layer, index);
3893 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3894 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3895 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3896 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3900 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
3903 // Take into account the mis-alignment (bring track to cluster plane)
3904 Double_t xTrOrig=track->GetX();
3905 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3906 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3907 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3908 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3910 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3913 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3914 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3915 c.SetSigmaYZ(track->GetSigmaYZ(layer));
3918 Int_t updated = track->UpdateMI(&c,chi2,index);
3919 // Bring the track back to detector plane in ideal geometry
3920 if (!track->Propagate(xTrOrig)) return 0;
3922 if(!updated) AliDebug(2,"update failed");
3926 //------------------------------------------------------------------------
3927 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3930 //DCA sigmas parameterization
3931 //to be paramterized using external parameters in future
3934 Double_t curv=track->GetC();
3935 sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
3936 sigmaz = 0.0110+4.37*TMath::Abs(curv);
3938 //------------------------------------------------------------------------
3939 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3942 // Clusters from delta electrons?
3944 Int_t entries = clusterArray->GetEntriesFast();
3945 if (entries<4) return;
3946 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3947 Int_t layer = cluster->GetLayer();
3948 if (layer>1) return;
3950 Int_t ncandidates=0;
3951 Float_t r = (layer>0)? 7:4;
3953 for (Int_t i=0;i<entries;i++){
3954 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3955 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3956 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3957 index[ncandidates] = i; //candidate to belong to delta electron track
3959 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3960 cl0->SetDeltaProbability(1);
3966 for (Int_t i=0;i<ncandidates;i++){
3967 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3968 if (cl0->GetDeltaProbability()>0.8) continue;
3971 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3972 sumy=sumz=sumy2=sumyz=sumw=0.0;
3973 for (Int_t j=0;j<ncandidates;j++){
3975 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3977 Float_t dz = cl0->GetZ()-cl1->GetZ();
3978 Float_t dy = cl0->GetY()-cl1->GetY();
3979 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3980 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3981 y[ncl] = cl1->GetY();
3982 z[ncl] = cl1->GetZ();
3983 sumy+= y[ncl]*weight;
3984 sumz+= z[ncl]*weight;
3985 sumy2+=y[ncl]*y[ncl]*weight;
3986 sumyz+=y[ncl]*z[ncl]*weight;
3991 if (ncl<4) continue;
3992 Float_t det = sumw*sumy2 - sumy*sumy;
3994 if (TMath::Abs(det)>0.01){
3995 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3996 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3997 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4000 Float_t z0 = sumyz/sumy;
4001 delta = TMath::Abs(cl0->GetZ()-z0);
4004 cl0->SetDeltaProbability(1-20.*delta);
4008 //------------------------------------------------------------------------
4009 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4014 track->UpdateESDtrack(flags);
4015 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4016 if (oldtrack) delete oldtrack;
4017 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4018 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4019 // printf("Problem\n");
4022 //------------------------------------------------------------------------
4023 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4025 // Get nearest upper layer close to the point xr.
4026 // rough approximation
4028 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4029 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4031 for (Int_t i=0;i<6;i++){
4032 if (radius<kRadiuses[i]){
4039 //------------------------------------------------------------------------
4040 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4041 //--------------------------------------------------------------------
4042 // Fill a look-up table with mean material
4043 //--------------------------------------------------------------------
4047 Double_t point1[3],point2[3];
4048 Double_t phi,cosphi,sinphi,z;
4049 // 0-5 layers, 6 pipe, 7-8 shields
4050 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4051 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4053 Int_t ifirst=0,ilast=0;
4054 if(material.Contains("Pipe")) {
4056 } else if(material.Contains("Shields")) {
4058 } else if(material.Contains("Layers")) {
4061 Error("BuildMaterialLUT","Wrong layer name\n");
4064 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4065 Double_t param[5]={0.,0.,0.,0.,0.};
4066 for (Int_t i=0; i<n; i++) {
4067 phi = 2.*TMath::Pi()*gRandom->Rndm();
4068 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4069 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4070 point1[0] = rmin[imat]*cosphi;
4071 point1[1] = rmin[imat]*sinphi;
4073 point2[0] = rmax[imat]*cosphi;
4074 point2[1] = rmax[imat]*sinphi;
4076 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4077 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4079 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4081 fxOverX0Layer[imat] = param[1];
4082 fxTimesRhoLayer[imat] = param[0]*param[4];
4083 } else if(imat==6) {
4084 fxOverX0Pipe = param[1];
4085 fxTimesRhoPipe = param[0]*param[4];
4086 } else if(imat==7) {
4087 fxOverX0Shield[0] = param[1];
4088 fxTimesRhoShield[0] = param[0]*param[4];
4089 } else if(imat==8) {
4090 fxOverX0Shield[1] = param[1];
4091 fxTimesRhoShield[1] = param[0]*param[4];
4095 printf("%s\n",material.Data());
4096 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4097 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4098 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4099 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4100 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4101 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4102 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4103 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4104 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4108 //------------------------------------------------------------------------
4109 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4110 TString direction) {
4111 //-------------------------------------------------------------------
4112 // Propagate beyond beam pipe and correct for material
4113 // (material budget in different ways according to fUseTGeo value)
4114 // Add time if going outward (PropagateTo or PropagateToTGeo)
4115 //-------------------------------------------------------------------
4117 // Define budget mode:
4118 // 0: material from AliITSRecoParam (hard coded)
4119 // 1: material from TGeo in one step (on the fly)
4120 // 2: material from lut
4121 // 3: material from TGeo in one step (same for all hypotheses)
4134 if(fTrackingPhase.Contains("Clusters2Tracks"))
4135 { mode=3; } else { mode=1; }
4138 if(fTrackingPhase.Contains("Clusters2Tracks"))
4139 { mode=3; } else { mode=2; }
4145 if(fTrackingPhase.Contains("Default")) mode=0;
4147 Int_t index=fCurrentEsdTrack;
4149 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4150 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4152 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4154 Double_t xOverX0,x0,lengthTimesMeanDensity;
4158 xOverX0 = AliITSRecoParam::GetdPipe();
4159 x0 = AliITSRecoParam::GetX0Be();
4160 lengthTimesMeanDensity = xOverX0*x0;
4161 lengthTimesMeanDensity *= dir;
4162 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4165 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4168 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4169 xOverX0 = fxOverX0Pipe;
4170 lengthTimesMeanDensity = fxTimesRhoPipe;
4171 lengthTimesMeanDensity *= dir;
4172 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4175 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4176 if(fxOverX0PipeTrks[index]<0) {
4177 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4178 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4179 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4180 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4181 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4184 xOverX0 = fxOverX0PipeTrks[index];
4185 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4186 lengthTimesMeanDensity *= dir;
4187 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4193 //------------------------------------------------------------------------
4194 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4196 TString direction) {
4197 //-------------------------------------------------------------------
4198 // Propagate beyond SPD or SDD shield and correct for material
4199 // (material budget in different ways according to fUseTGeo value)
4200 // Add time if going outward (PropagateTo or PropagateToTGeo)
4201 //-------------------------------------------------------------------
4203 // Define budget mode:
4204 // 0: material from AliITSRecoParam (hard coded)
4205 // 1: material from TGeo in steps of X cm (on the fly)
4206 // X = AliITSRecoParam::GetStepSizeTGeo()
4207 // 2: material from lut
4208 // 3: material from TGeo in one step (same for all hypotheses)
4221 if(fTrackingPhase.Contains("Clusters2Tracks"))
4222 { mode=3; } else { mode=1; }
4225 if(fTrackingPhase.Contains("Clusters2Tracks"))
4226 { mode=3; } else { mode=2; }
4232 if(fTrackingPhase.Contains("Default")) mode=0;
4234 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4236 Int_t shieldindex=0;
4237 if (shield.Contains("SDD")) { // SDDouter
4238 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4240 } else if (shield.Contains("SPD")) { // SPDouter
4241 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4244 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4248 // do nothing if we are already beyond the shield
4249 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4250 if(dir<0 && rTrack > rToGo) return 1; // going outward
4251 if(dir>0 && rTrack < rToGo) return 1; // going inward
4255 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4257 Int_t index=2*fCurrentEsdTrack+shieldindex;
4259 Double_t xOverX0,x0,lengthTimesMeanDensity;
4264 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4265 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4266 lengthTimesMeanDensity = xOverX0*x0;
4267 lengthTimesMeanDensity *= dir;
4268 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4271 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4272 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4275 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4276 xOverX0 = fxOverX0Shield[shieldindex];
4277 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4278 lengthTimesMeanDensity *= dir;
4279 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4282 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4283 if(fxOverX0ShieldTrks[index]<0) {
4284 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4285 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4286 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4287 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4288 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4291 xOverX0 = fxOverX0ShieldTrks[index];
4292 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4293 lengthTimesMeanDensity *= dir;
4294 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4300 //------------------------------------------------------------------------
4301 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4303 Double_t oldGlobXYZ[3],
4304 TString direction) {
4305 //-------------------------------------------------------------------
4306 // Propagate beyond layer and correct for material
4307 // (material budget in different ways according to fUseTGeo value)
4308 // Add time if going outward (PropagateTo or PropagateToTGeo)
4309 //-------------------------------------------------------------------
4311 // Define budget mode:
4312 // 0: material from AliITSRecoParam (hard coded)
4313 // 1: material from TGeo in stepsof X cm (on the fly)
4314 // X = AliITSRecoParam::GetStepSizeTGeo()
4315 // 2: material from lut
4316 // 3: material from TGeo in one step (same for all hypotheses)
4329 if(fTrackingPhase.Contains("Clusters2Tracks"))
4330 { mode=3; } else { mode=1; }
4333 if(fTrackingPhase.Contains("Clusters2Tracks"))
4334 { mode=3; } else { mode=2; }
4340 if(fTrackingPhase.Contains("Default")) mode=0;
4342 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4344 Double_t r=fgLayers[layerindex].GetR();
4345 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4347 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4349 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4351 Int_t index=6*fCurrentEsdTrack+layerindex;
4354 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4357 // back before material (no correction)
4359 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4360 if (!t->GetLocalXat(rOld,xOld)) return 0;
4361 if (!t->Propagate(xOld)) return 0;
4365 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4366 lengthTimesMeanDensity = xOverX0*x0;
4367 lengthTimesMeanDensity *= dir;
4368 // Bring the track beyond the material
4369 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4372 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4373 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4376 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4377 xOverX0 = fxOverX0Layer[layerindex];
4378 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4379 lengthTimesMeanDensity *= dir;
4380 // Bring the track beyond the material
4381 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4384 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4385 if(fxOverX0LayerTrks[index]<0) {
4386 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4387 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4388 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4389 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4390 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4391 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4394 xOverX0 = fxOverX0LayerTrks[index];
4395 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4396 lengthTimesMeanDensity *= dir;
4397 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4404 //------------------------------------------------------------------------
4405 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4406 //-----------------------------------------------------------------
4407 // Initialize LUT for storing material for each prolonged track
4408 //-----------------------------------------------------------------
4409 fxOverX0PipeTrks = new Float_t[ntracks];
4410 fxTimesRhoPipeTrks = new Float_t[ntracks];
4411 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4412 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4413 fxOverX0LayerTrks = new Float_t[ntracks*6];
4414 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4416 for(Int_t i=0; i<ntracks; i++) {
4417 fxOverX0PipeTrks[i] = -1.;
4418 fxTimesRhoPipeTrks[i] = -1.;
4420 for(Int_t j=0; j<ntracks*2; j++) {
4421 fxOverX0ShieldTrks[j] = -1.;
4422 fxTimesRhoShieldTrks[j] = -1.;
4424 for(Int_t k=0; k<ntracks*6; k++) {
4425 fxOverX0LayerTrks[k] = -1.;
4426 fxTimesRhoLayerTrks[k] = -1.;
4433 //------------------------------------------------------------------------
4434 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4435 //-----------------------------------------------------------------
4436 // Delete LUT for storing material for each prolonged track
4437 //-----------------------------------------------------------------
4438 if(fxOverX0PipeTrks) {
4439 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4441 if(fxOverX0ShieldTrks) {
4442 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4445 if(fxOverX0LayerTrks) {
4446 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4448 if(fxTimesRhoPipeTrks) {
4449 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4451 if(fxTimesRhoShieldTrks) {
4452 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4454 if(fxTimesRhoLayerTrks) {
4455 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4459 //------------------------------------------------------------------------
4460 void AliITStrackerMI::SetForceSkippingOfLayer() {
4461 //-----------------------------------------------------------------
4462 // Check if we are forced to skip layers
4463 // either we set to skip them in RecoParam
4464 // or they were off during data-taking
4465 //-----------------------------------------------------------------
4467 const AliEventInfo *eventInfo = GetEventInfo();
4469 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4470 fForceSkippingOfLayer[l] = 0;
4472 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4476 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4477 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4479 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4480 } else if(l==2 || l==3) {
4481 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4483 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4489 //------------------------------------------------------------------------
4490 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4491 Int_t ilayer,Int_t idet) const {
4492 //-----------------------------------------------------------------
4493 // This method is used to decide whether to allow a prolongation
4494 // without clusters, because we want to skip the layer.
4495 // In this case the return value is > 0:
4496 // return 1: the user requested to skip a layer
4497 // return 2: track outside z acceptance
4498 //-----------------------------------------------------------------
4500 if (ForceSkippingOfLayer(ilayer)) return 1;
4502 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4504 if (idet<0 && // out in z
4505 ilayer>innerLayCanSkip &&
4506 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4507 // check if track will cross SPD outer layer
4508 Double_t phiAtSPD2,zAtSPD2;
4509 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4510 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4512 return 2; // always allow skipping, changed on 05.11.2009
4517 //------------------------------------------------------------------------
4518 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4519 Int_t ilayer,Int_t idet,
4520 Double_t dz,Double_t dy,
4521 Bool_t noClusters) const {
4522 //-----------------------------------------------------------------
4523 // This method is used to decide whether to allow a prolongation
4524 // without clusters, because there is a dead zone in the road.
4525 // In this case the return value is > 0:
4526 // return 1: dead zone at z=0,+-7cm in SPD
4527 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4528 // return 2: all road is "bad" (dead or noisy) from the OCDB
4529 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4530 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4531 //-----------------------------------------------------------------
4533 // check dead zones at z=0,+-7cm in the SPD
4534 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4535 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4536 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4537 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4538 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4539 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4540 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4541 for (Int_t i=0; i<3; i++)
4542 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4543 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4544 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4548 // check bad zones from OCDB
4549 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4551 if (idet<0) return 0;
4553 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4556 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4557 if (ilayer==0 || ilayer==1) { // ---------- SPD
4559 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4561 detSizeFactorX *= 2.;
4562 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4565 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4566 if (detType==2) segm->SetLayer(ilayer+1);
4567 Float_t detSizeX = detSizeFactorX*segm->Dx();
4568 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4570 // check if the road overlaps with bad chips
4572 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4573 Float_t zlocmin = zloc-dz;
4574 Float_t zlocmax = zloc+dz;
4575 Float_t xlocmin = xloc-dy;
4576 Float_t xlocmax = xloc+dy;
4577 Int_t chipsInRoad[100];
4579 // check if road goes out of detector
4580 Bool_t touchNeighbourDet=kFALSE;
4581 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4582 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4583 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4584 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4585 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));
4587 // check if this detector is bad
4589 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4590 if(!touchNeighbourDet) {
4591 return 2; // all detectors in road are bad
4593 return 3; // at least one is bad
4597 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4598 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4599 if (!nChipsInRoad) return 0;
4601 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4602 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4603 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4604 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4605 if (det.IsChipBad(chipsInRoad[iCh])) {
4613 if(!touchNeighbourDet) {
4614 AliDebug(2,"all bad in road");
4615 return 2; // all chips in road are bad
4617 return 3; // at least a bad chip in road
4622 AliDebug(2,"at least a bad in road");
4623 return 3; // at least a bad chip in road
4627 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4628 || !noClusters) return 0;
4630 // There are no clusters in road: check if there is at least
4631 // a bad SPD pixel or SDD anode or SSD strips on both sides
4633 Int_t idetInITS=idet;
4634 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4636 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4637 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4640 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4644 //------------------------------------------------------------------------
4645 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4646 const AliITStrackMI *track,
4647 Float_t &xloc,Float_t &zloc) const {
4648 //-----------------------------------------------------------------
4649 // Gives position of track in local module ref. frame
4650 //-----------------------------------------------------------------
4655 if(idet<0) return kTRUE; // track out of z acceptance of layer
4657 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4659 Int_t lad = Int_t(idet/ndet) + 1;
4661 Int_t det = idet - (lad-1)*ndet + 1;
4663 Double_t xyzGlob[3],xyzLoc[3];
4665 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4666 // take into account the misalignment: xyz at real detector plane
4667 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4669 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4671 xloc = (Float_t)xyzLoc[0];
4672 zloc = (Float_t)xyzLoc[2];
4676 //------------------------------------------------------------------------
4677 //------------------------------------------------------------------------
4678 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4680 // Method to be optimized further:
4681 // Aim: decide whether a track can be used for PlaneEff evaluation
4682 // the decision is taken based on the track quality at the layer under study
4683 // no information on the clusters on this layer has to be used
4684 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4685 // the cut is done on number of sigmas from the boundaries
4687 // Input: Actual track, layer [0,5] under study
4689 // Return: kTRUE if this is a good track
4691 // it will apply a pre-selection to obtain good quality tracks.
4692 // Here also you will have the possibility to put a control on the
4693 // impact point of the track on the basic block, in order to exclude border regions
4694 // this will be done by calling a proper method of the AliITSPlaneEff class.
4696 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4697 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4699 Int_t index[AliITSgeomTGeo::kNLayers];
4701 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4703 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4704 index[k]=clusters[k];
4708 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4709 AliITSlayer &layer=fgLayers[ilayer];
4710 Double_t r=layer.GetR();
4711 AliITStrackMI tmp(*track);
4713 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4714 Int_t ncl_out=0; Int_t ncl_in=0;
4715 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4716 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4717 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4718 // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4719 if(index[lay]>=0)ncl_out++;
4721 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4722 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4723 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4724 if (index[lay]>=0) ncl_in++;
4726 Int_t ncl=ncl_out+ncl_in;
4727 Bool_t nextout = kFALSE;
4728 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4729 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4730 Bool_t nextin = kFALSE;
4731 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4732 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4733 // maximum number of missing clusters allowed in outermost layers
4734 if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4736 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4737 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4739 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4740 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4741 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4742 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4746 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4747 Int_t idet=layer.FindDetectorIndex(phi,z);
4748 if(idet<0) { AliInfo(Form("cannot find detector"));
4751 // here check if it has good Chi Square.
4753 //propagate to the intersection with the detector plane
4754 const AliITSdetector &det=layer.GetDetector(idet);
4755 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4759 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4760 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4761 if(key>fPlaneEff->Nblock()) return kFALSE;
4762 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4763 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4765 // DEFINITION OF SEARCH ROAD FOR accepting a track
4767 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4768 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4769 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4770 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4771 // done for RecPoints
4773 // exclude tracks at boundary between detectors
4774 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4775 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4776 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4777 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4778 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4779 if ( (locx-dx < blockXmn+boundaryWidth) ||
4780 (locx+dx > blockXmx-boundaryWidth) ||
4781 (locz-dz < blockZmn+boundaryWidth) ||
4782 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4785 //------------------------------------------------------------------------
4786 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4788 // This Method has to be optimized! For the time-being it uses the same criteria
4789 // as those used in the search of extra clusters for overlapping modules.
4791 // Method Purpose: estabilish whether a track has produced a recpoint or not
4792 // in the layer under study (For Plane efficiency)
4794 // inputs: AliITStrackMI* track (pointer to a usable track)
4796 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4797 // traversed by this very track. In details:
4798 // - if a cluster can be associated to the track then call
4799 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4800 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4803 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4804 AliITSlayer &layer=fgLayers[ilayer];
4805 Double_t r=layer.GetR();
4806 AliITStrackMI tmp(*track);
4810 if (!tmp.GetPhiZat(r,phi,z)) return;
4811 Int_t idet=layer.FindDetectorIndex(phi,z);
4813 if(idet<0) { AliInfo(Form("cannot find detector"));
4817 //propagate to the intersection with the detector plane
4818 const AliITSdetector &det=layer.GetDetector(idet);
4819 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4823 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4825 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4826 TMath::Sqrt(tmp.GetSigmaZ2() +
4827 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4828 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4829 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4830 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4831 TMath::Sqrt(tmp.GetSigmaY2() +
4832 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4833 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4834 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4836 // road in global (rphi,z) [i.e. in tracking ref. system]
4837 Double_t zmin = tmp.GetZ() - dz;
4838 Double_t zmax = tmp.GetZ() + dz;
4839 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4840 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4842 // select clusters in road
4843 layer.SelectClusters(zmin,zmax,ymin,ymax);
4845 // Define criteria for track-cluster association
4846 Double_t msz = tmp.GetSigmaZ2() +
4847 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4848 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4849 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4850 Double_t msy = tmp.GetSigmaY2() +
4851 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4852 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4853 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4854 if (tmp.GetConstrain()) {
4855 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4856 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4858 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4859 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4861 msz = 1./msz; // 1/RoadZ^2
4862 msy = 1./msy; // 1/RoadY^2
4865 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4867 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4868 //Double_t tolerance=0.2;
4869 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4870 idetc = cl->GetDetectorIndex();
4871 if(idet!=idetc) continue;
4872 //Int_t ilay = cl->GetLayer();
4874 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4875 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4877 Double_t chi2=tmp.GetPredictedChi2(cl);
4878 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4882 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4884 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4885 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4886 if(key>fPlaneEff->Nblock()) return;
4887 Bool_t found=kFALSE;
4890 while ((cl=layer.GetNextCluster(clidx))!=0) {
4891 idetc = cl->GetDetectorIndex();
4892 if(idet!=idetc) continue;
4893 // here real control to see whether the cluster can be associated to the track.
4894 // cluster not associated to track
4895 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4896 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
4897 // calculate track-clusters chi2
4898 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4899 // in particular, the error associated to the cluster
4900 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4902 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4904 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4905 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4906 // track->SetExtraModule(ilayer,idetExtra);
4908 if(!fPlaneEff->UpDatePlaneEff(found,key))
4909 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4910 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4911 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
4912 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4913 Int_t cltype[2]={-999,-999};
4916 Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
4920 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4921 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4924 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4925 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4926 cltype[0]=layer.GetCluster(ci)->GetNy();
4927 cltype[1]=layer.GetCluster(ci)->GetNz();
4929 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4930 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
4931 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4932 // It is computed properly by calling the method
4933 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4935 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4936 //if (tmp.PropagateTo(x,0.,0.)) {
4937 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4938 AliCluster c(*layer.GetCluster(ci));
4939 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4940 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4941 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
4942 clu[2]=TMath::Sqrt(c.GetSigmaY2());
4943 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4946 // Compute the angles between the track and the module
4947 // compute the angle "in phi direction", i.e. the angle in the transverse plane
4948 // between the normal to the module and the projection (in the transverse plane) of the
4950 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
4951 Float_t tgl = tmp.GetTgl();
4952 Float_t phitr = tmp.GetSnp();
4953 phitr = TMath::ASin(phitr);
4954 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
4956 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
4957 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
4959 alpha = tmp.GetAlpha();
4960 Double_t phiglob = alpha+phitr;
4962 p[0] = TMath::Cos(phiglob);
4963 p[1] = TMath::Sin(phiglob);
4965 TVector3 pvec(p[0],p[1],p[2]);
4966 TVector3 normvec(rot[1],rot[4],rot[7]);
4967 Double_t angle = pvec.Angle(normvec);
4969 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
4970 angle *= 180./TMath::Pi();
4973 TVector3 pt(p[0],p[1],0);
4974 TVector3 normt(rot[1],rot[4],0);
4975 Double_t anglet = pt.Angle(normt);
4977 Double_t phiPt = TMath::ATan2(p[1],p[0]);
4978 if(phiPt<0)phiPt+=2.*TMath::Pi();
4979 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
4980 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
4981 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
4982 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
4983 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
4984 anglet *= 180./TMath::Pi();
4986 AngleModTrack[2]=(Float_t) angle;
4987 AngleModTrack[0]=(Float_t) anglet;
4988 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
4989 AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
4990 AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
4991 AngleModTrack[1]*=180./TMath::Pi(); // in degree
4993 fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);