1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------------------------
18 // Implementation of the ITS tracker class
19 // It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
20 // and fills with them the ESD
21 // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
22 // Current support and development:
23 // Andrea Dainese, andrea.dainese@lnl.infn.it
24 // dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25 // Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26 // Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
27 //-------------------------------------------------------------------------
31 #include <TTreeStream.h>
32 #include <TDatabasePDG.h>
37 #include "AliESDEvent.h"
38 #include "AliESDtrack.h"
39 #include "AliESDVertex.h"
42 #include "AliITSRecPoint.h"
43 #include "AliITSgeomTGeo.h"
44 #include "AliITSReconstructor.h"
45 #include "AliTrackPointArray.h"
46 #include "AliAlignObj.h"
47 #include "AliITSClusterParam.h"
48 #include "AliCDBManager.h"
49 #include "AliCDBEntry.h"
50 #include "AliITSsegmentation.h"
51 #include "AliITSCalibration.h"
52 #include "AliITSCalibrationSPD.h"
53 #include "AliITSCalibrationSDD.h"
54 #include "AliITSCalibrationSSD.h"
55 #include "AliITSPlaneEff.h"
56 #include "AliITSPlaneEffSPD.h"
57 #include "AliITSPlaneEffSDD.h"
58 #include "AliITSPlaneEffSSD.h"
59 #include "AliITStrackerMI.h"
61 ClassImp(AliITStrackerMI)
63 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
65 AliITStrackerMI::AliITStrackerMI():AliTracker(),
75 fLastLayerToTrackTo(0),
78 fTrackingPhase("Default"),
84 fxTimesRhoPipeTrks(0),
85 fxOverX0ShieldTrks(0),
86 fxTimesRhoShieldTrks(0),
88 fxTimesRhoLayerTrks(0),
95 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
96 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
97 for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
99 //------------------------------------------------------------------------
100 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
101 fI(AliITSgeomTGeo::GetNLayers()),
110 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
113 fTrackingPhase("Default"),
119 fxTimesRhoPipeTrks(0),
120 fxOverX0ShieldTrks(0),
121 fxTimesRhoShieldTrks(0),
122 fxOverX0LayerTrks(0),
123 fxTimesRhoLayerTrks(0),
125 fITSChannelStatus(0),
128 //--------------------------------------------------------------------
129 //This is the AliITStrackerMI constructor
130 //--------------------------------------------------------------------
132 AliWarning("\"geom\" is actually a dummy argument !");
138 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
139 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
140 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
142 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
143 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
144 Double_t poff=TMath::ATan2(y,x);
146 Double_t r=TMath::Sqrt(x*x + y*y);
148 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
149 r += TMath::Sqrt(x*x + y*y);
150 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
151 r += TMath::Sqrt(x*x + y*y);
152 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
153 r += TMath::Sqrt(x*x + y*y);
156 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
158 for (Int_t j=1; j<nlad+1; j++) {
159 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
160 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
161 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
163 Double_t txyz[3]={0.};
164 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
165 m.LocalToMaster(txyz,xyz);
166 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
167 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
169 if (phi<0) phi+=TMath::TwoPi();
170 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
172 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
173 new(&det) AliITSdetector(r,phi);
174 // compute the real radius (with misalignment)
175 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
177 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
178 mmisal.LocalToMaster(txyz,xyz);
179 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
180 det.SetRmisal(rmisal);
182 } // end loop on detectors
183 } // end loop on ladders
184 } // end loop on layers
187 fI=AliITSgeomTGeo::GetNLayers();
190 fConstraint[0]=1; fConstraint[1]=0;
192 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
193 AliITSReconstructor::GetRecoParam()->GetYVdef(),
194 AliITSReconstructor::GetRecoParam()->GetZVdef()};
195 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
196 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
197 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
198 SetVertex(xyzVtx,ersVtx);
200 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=AliITSRecoParam::GetLayersNotToSkip(i);
201 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
202 for (Int_t i=0;i<100000;i++){
203 fBestTrackIndex[i]=0;
206 // store positions of centre of SPD modules (in z)
208 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
209 fSPDdetzcentre[0] = tr[2];
210 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
211 fSPDdetzcentre[1] = tr[2];
212 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
213 fSPDdetzcentre[2] = tr[2];
214 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
215 fSPDdetzcentre[3] = tr[2];
217 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
218 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
219 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
223 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
224 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
226 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
228 // only for plane efficiency evaluation
229 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff()) {
230 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
231 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane))
232 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
233 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
234 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
235 else fPlaneEff = new AliITSPlaneEffSSD();
236 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
237 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
238 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
241 //------------------------------------------------------------------------
242 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
244 fBestTrack(tracker.fBestTrack),
245 fTrackToFollow(tracker.fTrackToFollow),
246 fTrackHypothesys(tracker.fTrackHypothesys),
247 fBestHypothesys(tracker.fBestHypothesys),
248 fOriginal(tracker.fOriginal),
249 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
250 fPass(tracker.fPass),
251 fAfterV0(tracker.fAfterV0),
252 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
253 fCoefficients(tracker.fCoefficients),
255 fTrackingPhase(tracker.fTrackingPhase),
256 fUseTGeo(tracker.fUseTGeo),
257 fNtracks(tracker.fNtracks),
258 fxOverX0Pipe(tracker.fxOverX0Pipe),
259 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
261 fxTimesRhoPipeTrks(0),
262 fxOverX0ShieldTrks(0),
263 fxTimesRhoShieldTrks(0),
264 fxOverX0LayerTrks(0),
265 fxTimesRhoLayerTrks(0),
266 fDebugStreamer(tracker.fDebugStreamer),
267 fITSChannelStatus(tracker.fITSChannelStatus),
268 fDetTypeRec(tracker.fDetTypeRec),
269 fPlaneEff(tracker.fPlaneEff) {
273 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
276 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
277 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
280 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
281 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
284 //------------------------------------------------------------------------
285 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
286 //Assignment operator
287 this->~AliITStrackerMI();
288 new(this) AliITStrackerMI(tracker);
291 //------------------------------------------------------------------------
292 AliITStrackerMI::~AliITStrackerMI()
297 if (fCoefficients) delete [] fCoefficients;
298 DeleteTrksMaterialLUT();
299 if (fDebugStreamer) {
300 //fDebugStreamer->Close();
301 delete fDebugStreamer;
303 if(fITSChannelStatus) delete fITSChannelStatus;
304 if(fPlaneEff) delete fPlaneEff;
306 //------------------------------------------------------------------------
307 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
308 //--------------------------------------------------------------------
309 //This function set masks of the layers which must be not skipped
310 //--------------------------------------------------------------------
311 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fLayersNotToSkip[i]=l[i];
313 //------------------------------------------------------------------------
314 void AliITStrackerMI::ReadBadFromDetTypeRec() {
315 //--------------------------------------------------------------------
316 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
318 //--------------------------------------------------------------------
320 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
322 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels\n");
324 if(!fDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
327 if(fITSChannelStatus) delete fITSChannelStatus;
328 fITSChannelStatus = new AliITSChannelStatus(AliCDBManager::Instance());
330 // ITS detectors and chips
331 Int_t i=0,j=0,k=0,ndet=0;
332 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
333 ndet=AliITSgeomTGeo::GetNDetectors(i);
334 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
335 for (k=1; k<ndet+1; k++) {
336 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
337 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fDetTypeRec);
338 } // end loop on detectors
339 } // end loop on ladders
340 } // end loop on layers
344 //------------------------------------------------------------------------
345 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
346 //--------------------------------------------------------------------
347 //This function loads ITS clusters
348 //--------------------------------------------------------------------
349 TBranch *branch=cTree->GetBranch("ITSRecPoints");
351 Error("LoadClusters"," can't get the branch !\n");
355 static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
356 branch->SetAddress(&clusters);
358 Int_t i=0,j=0,ndet=0;
360 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
361 ndet=fgLayers[i].GetNdetectors();
362 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
363 for (; j<jmax; j++) {
364 if (!cTree->GetEvent(j)) continue;
365 Int_t ncl=clusters->GetEntriesFast();
366 SignDeltas(clusters,GetZ());
369 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
370 detector=c->GetDetectorIndex();
372 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
374 fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
377 // add dead zone "virtual" cluster in SPD, if there is a cluster within
378 // zwindow cm from the dead zone
379 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
380 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
381 Int_t lab[4] = {0,0,0,detector};
382 Int_t info[3] = {0,0,i};
383 Float_t q = 0.; // this identifies virtual clusters
384 Float_t hit[5] = {xdead,
386 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
387 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
389 Bool_t local = kTRUE;
390 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
391 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
392 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
393 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
394 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
395 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
396 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
397 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
398 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
399 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
400 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
401 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
402 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
403 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
404 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
405 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
406 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
407 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
408 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
410 } // "virtual" clusters in SPD
414 fgLayers[i].ResetRoad(); //road defined by the cluster density
415 fgLayers[i].SortClusters();
422 //------------------------------------------------------------------------
423 void AliITStrackerMI::UnloadClusters() {
424 //--------------------------------------------------------------------
425 //This function unloads ITS clusters
426 //--------------------------------------------------------------------
427 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
429 //------------------------------------------------------------------------
430 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
431 //--------------------------------------------------------------------
432 // Publishes all pointers to clusters known to the tracker into the
433 // passed object array.
434 // The ownership is not transfered - the caller is not expected to delete
436 //--------------------------------------------------------------------
438 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
439 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
440 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
447 //------------------------------------------------------------------------
448 static Int_t CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
449 //--------------------------------------------------------------------
450 // Correction for the material between the TPC and the ITS
451 //--------------------------------------------------------------------
452 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
453 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
454 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
455 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
456 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
457 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
458 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
459 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
461 Error("CorrectForTPCtoITSDeadZoneMaterial","Track is already in the dead zone !");
467 //------------------------------------------------------------------------
468 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
469 //--------------------------------------------------------------------
470 // This functions reconstructs ITS tracks
471 // The clusters must be already loaded !
472 //--------------------------------------------------------------------
475 fTrackingPhase="Clusters2Tracks";
477 TObjArray itsTracks(15000);
479 fEsd = event; // store pointer to the esd
481 // temporary (for cosmics)
482 if(event->GetVertex()) {
483 TString title = event->GetVertex()->GetTitle();
484 if(title.Contains("cosmics")) {
485 Double_t xyz[3]={GetX(),GetY(),GetZ()};
486 Double_t exyz[3]={0.1,0.1,0.1};
492 {/* Read ESD tracks */
493 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
494 Int_t nentr=event->GetNumberOfTracks();
495 Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
497 AliESDtrack *esd=event->GetTrack(nentr);
499 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
500 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
501 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
502 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
505 t=new AliITStrackMI(*esd);
506 } catch (const Char_t *msg) {
507 //Warning("Clusters2Tracks",msg);
511 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
512 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
515 // look at the ESD mass hypothesys !
516 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
518 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
520 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
521 //track - can be V0 according to TPC
523 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
527 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
531 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
535 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
540 t->SetReconstructed(kFALSE);
541 itsTracks.AddLast(t);
542 fOriginal.AddLast(t);
544 } /* End Read ESD tracks */
548 Int_t nentr=itsTracks.GetEntriesFast();
549 fTrackHypothesys.Expand(nentr);
550 fBestHypothesys.Expand(nentr);
551 MakeCoefficients(nentr);
552 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
554 // THE TWO TRACKING PASSES
555 for (fPass=0; fPass<2; fPass++) {
556 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
557 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
558 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
559 if (t==0) continue; //this track has been already tracked
560 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
561 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
562 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
563 if (fConstraint[fPass]) {
564 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
565 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
568 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
570 ResetTrackToFollow(*t);
573 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
576 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
578 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
579 if (!besttrack) continue;
580 besttrack->SetLabel(tpcLabel);
581 // besttrack->CookdEdx();
583 besttrack->SetFakeRatio(1.);
584 CookLabel(besttrack,0.); //For comparison only
585 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
587 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
589 t->SetReconstructed(kTRUE);
592 GetBestHypothesysMIP(itsTracks);
593 } // end loop on the two tracking passes
595 if(event->GetNumberOfV0s()>0) UpdateTPCV0(event);
596 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) FindV02(event);
601 Int_t entries = fTrackHypothesys.GetEntriesFast();
602 for (Int_t ientry=0; ientry<entries; ientry++) {
603 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
604 if (array) array->Delete();
605 delete fTrackHypothesys.RemoveAt(ientry);
608 fTrackHypothesys.Delete();
609 fBestHypothesys.Delete();
611 delete [] fCoefficients;
613 DeleteTrksMaterialLUT();
615 Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
617 fTrackingPhase="Default";
621 //------------------------------------------------------------------------
622 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
623 //--------------------------------------------------------------------
624 // This functions propagates reconstructed ITS tracks back
625 // The clusters must be loaded !
626 //--------------------------------------------------------------------
627 fTrackingPhase="PropagateBack";
628 Int_t nentr=event->GetNumberOfTracks();
629 Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
632 for (Int_t i=0; i<nentr; i++) {
633 AliESDtrack *esd=event->GetTrack(i);
635 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
636 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
640 t=new AliITStrackMI(*esd);
641 } catch (const Char_t *msg) {
642 //Warning("PropagateBack",msg);
646 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
648 ResetTrackToFollow(*t);
650 // propagate to vertex [SR, GSI 17.02.2003]
651 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
652 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
653 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
654 fTrackToFollow.StartTimeIntegral();
655 // from vertex to outside pipe
656 CorrectForPipeMaterial(&fTrackToFollow,"outward");
659 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
660 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
661 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
665 fTrackToFollow.SetLabel(t->GetLabel());
666 //fTrackToFollow.CookdEdx();
667 CookLabel(&fTrackToFollow,0.); //For comparison only
668 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
669 //UseClusters(&fTrackToFollow);
675 Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
677 fTrackingPhase="Default";
681 //------------------------------------------------------------------------
682 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
683 //--------------------------------------------------------------------
684 // This functions refits ITS tracks using the
685 // "inward propagated" TPC tracks
686 // The clusters must be loaded !
687 //--------------------------------------------------------------------
688 fTrackingPhase="RefitInward";
689 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) RefitV02(event);
690 Int_t nentr=event->GetNumberOfTracks();
691 Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
694 for (Int_t i=0; i<nentr; i++) {
695 AliESDtrack *esd=event->GetTrack(i);
697 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
698 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
699 if (esd->GetStatus()&AliESDtrack::kTPCout)
700 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
704 t=new AliITStrackMI(*esd);
705 } catch (const Char_t *msg) {
706 //Warning("RefitInward",msg);
710 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
711 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
716 ResetTrackToFollow(*t);
717 fTrackToFollow.ResetClusters();
719 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
720 fTrackToFollow.ResetCovariance(10.);
723 Bool_t pe=AliITSReconstructor::GetRecoParam()->GetComputePlaneEff();
724 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
725 fTrackToFollow.SetLabel(t->GetLabel());
726 // fTrackToFollow.CookdEdx();
727 CookdEdx(&fTrackToFollow);
729 CookLabel(&fTrackToFollow,0.0); //For comparison only
732 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
733 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
734 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
735 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
736 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
737 Float_t r[3]={0.,0.,0.};
739 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
746 Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
748 fTrackingPhase="Default";
752 //------------------------------------------------------------------------
753 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
754 //--------------------------------------------------------------------
755 // Return pointer to a given cluster
756 //--------------------------------------------------------------------
757 Int_t l=(index & 0xf0000000) >> 28;
758 Int_t c=(index & 0x0fffffff) >> 00;
759 return fgLayers[l].GetCluster(c);
761 //------------------------------------------------------------------------
762 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
763 //--------------------------------------------------------------------
764 // Get track space point with index i
765 //--------------------------------------------------------------------
767 Int_t l=(index & 0xf0000000) >> 28;
768 Int_t c=(index & 0x0fffffff) >> 00;
769 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
770 Int_t idet = cl->GetDetectorIndex();
774 cl->GetGlobalXYZ(xyz);
775 cl->GetGlobalCov(cov);
777 p.SetCharge(cl->GetQ());
778 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
781 iLayer = AliGeomManager::kSPD1;
784 iLayer = AliGeomManager::kSPD2;
787 iLayer = AliGeomManager::kSDD1;
790 iLayer = AliGeomManager::kSDD2;
793 iLayer = AliGeomManager::kSSD1;
796 iLayer = AliGeomManager::kSSD2;
799 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
802 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
803 p.SetVolumeID((UShort_t)volid);
806 //------------------------------------------------------------------------
807 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
808 AliTrackPoint& p, const AliESDtrack *t) {
809 //--------------------------------------------------------------------
810 // Get track space point with index i
811 // (assign error estimated during the tracking)
812 //--------------------------------------------------------------------
814 Int_t l=(index & 0xf0000000) >> 28;
815 Int_t c=(index & 0x0fffffff) >> 00;
816 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
817 Int_t idet = cl->GetDetectorIndex();
818 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
820 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
822 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
823 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
824 Double_t alpha = t->GetAlpha();
825 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
826 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,AliTracker::GetBz()));
827 phi += alpha-det.GetPhi();
828 Float_t tgphi = TMath::Tan(phi);
830 Float_t tgl = t->GetTgl(); // tgl about const along track
831 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
833 Float_t errlocalx,errlocalz;
834 Bool_t addMisalErr=kFALSE;
835 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
839 cl->GetGlobalXYZ(xyz);
840 // cl->GetGlobalCov(cov);
841 Float_t pos[3] = {0.,0.,0.};
842 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
843 tmpcl.GetGlobalCov(cov);
846 p.SetCharge(cl->GetQ());
848 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
851 iLayer = AliGeomManager::kSPD1;
854 iLayer = AliGeomManager::kSPD2;
857 iLayer = AliGeomManager::kSDD1;
860 iLayer = AliGeomManager::kSDD2;
863 iLayer = AliGeomManager::kSSD1;
866 iLayer = AliGeomManager::kSSD2;
869 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
872 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
873 p.SetVolumeID((UShort_t)volid);
876 //------------------------------------------------------------------------
877 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
879 //--------------------------------------------------------------------
880 // Follow prolongation tree
881 //--------------------------------------------------------------------
883 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
884 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
887 AliESDtrack * esd = otrack->GetESDtrack();
888 if (esd->GetV0Index(0)>0) {
889 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
890 // mapping of ESD track is different as ITS track in Containers
891 // Need something more stable
892 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
893 for (Int_t i=0;i<3;i++){
894 Int_t index = esd->GetV0Index(i);
896 AliESDv0 * vertex = fEsd->GetV0(index);
897 if (vertex->GetStatus()<0) continue; // rejected V0
899 if (esd->GetSign()>0) {
900 vertex->SetIndex(0,esdindex);
902 vertex->SetIndex(1,esdindex);
906 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
908 bestarray = new TObjArray(5);
909 fBestHypothesys.AddAt(bestarray,esdindex);
913 //setup tree of the prolongations
915 static AliITStrackMI tracks[7][100];
916 AliITStrackMI *currenttrack;
917 static AliITStrackMI currenttrack1;
918 static AliITStrackMI currenttrack2;
919 static AliITStrackMI backuptrack;
921 Int_t nindexes[7][100];
922 Float_t normalizedchi2[100];
923 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
924 otrack->SetNSkipped(0);
925 new (&(tracks[6][0])) AliITStrackMI(*otrack);
927 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
928 Int_t modstatus = 1; // found
932 // follow prolongations
933 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
934 //printf("FollowProlongationTree: layer %d\n",ilayer);
937 AliITSlayer &layer=fgLayers[ilayer];
938 Double_t r = layer.GetR();
944 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
946 if (ntracks[ilayer]>=100) break;
947 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
948 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
949 if (ntracks[ilayer]>15+ilayer){
950 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
951 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
954 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
956 // material between SSD and SDD, SDD and SPD
958 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
960 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
964 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
965 Int_t idet=layer.FindDetectorIndex(phi,z);
967 Double_t trackGlobXYZ1[3];
968 currenttrack1.GetXYZ(trackGlobXYZ1);
970 // Get the budget to the primary vertex for the current track being prolonged
971 Double_t budgetToPrimVertex = GetEffectiveThickness();
973 // check if we allow a prolongation without point
974 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
976 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
977 // propagate to the layer radius
978 Double_t xToGo; vtrack->GetLocalXat(r,xToGo);
979 vtrack->AliExternalTrackParam::PropagateTo(xToGo,GetBz());
980 // apply correction for material of the current layer
981 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
982 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
983 vtrack->SetClIndex(ilayer,0);
984 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
985 LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc); // local module coords
986 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
987 if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
992 // track outside layer acceptance in z
993 if (idet<0) continue;
995 //propagate to the intersection with the detector plane
996 const AliITSdetector &det=layer.GetDetector(idet);
997 new(¤ttrack2) AliITStrackMI(currenttrack1);
998 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
999 currenttrack2.Propagate(det.GetPhi(),det.GetR());
1000 currenttrack1.SetDetectorIndex(idet);
1001 currenttrack2.SetDetectorIndex(idet);
1002 LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc); // local module coords
1005 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1007 // road in global (rphi,z) [i.e. in tracking ref. system]
1008 Double_t zmin,zmax,ymin,ymax;
1009 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1011 // select clusters in road
1012 layer.SelectClusters(zmin,zmax,ymin,ymax);
1013 //********************
1015 // Define criteria for track-cluster association
1016 Double_t msz = currenttrack1.GetSigmaZ2() +
1017 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1018 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1019 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1020 Double_t msy = currenttrack1.GetSigmaY2() +
1021 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1022 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1023 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1025 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1026 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1028 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1029 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1031 msz = 1./msz; // 1/RoadZ^2
1032 msy = 1./msy; // 1/RoadY^2
1036 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1038 const AliITSRecPoint *cl=0;
1040 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1041 Bool_t deadzoneSPD=kFALSE;
1042 currenttrack = ¤ttrack1;
1044 // check if the road contains a dead zone
1045 Bool_t noClusters = kFALSE;
1046 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1047 //if (noClusters) printf("no clusters in road\n");
1048 Double_t dz=0.5*(zmax-zmin);
1049 Double_t dy=0.5*(ymax-ymin);
1050 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1051 // create a prolongation without clusters (check also if there are no clusters in the road)
1054 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1055 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1056 updatetrack->SetClIndex(ilayer,0);
1058 modstatus = 5; // no cls in road
1059 } else if (dead==1) {
1060 modstatus = 7; // holes in z in SPD
1061 } else if (dead==2 || dead==3) {
1062 modstatus = 2; // dead from OCDB
1064 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1065 // apply correction for material of the current layer
1066 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1067 if (constrain) { // apply vertex constrain
1068 updatetrack->SetConstrain(constrain);
1069 Bool_t isPrim = kTRUE;
1070 if (ilayer<4) { // check that it's close to the vertex
1071 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1072 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1073 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1074 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1075 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1077 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1080 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1081 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1082 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1090 // loop over clusters in the road
1091 while ((cl=layer.GetNextCluster(clidx))!=0) {
1092 if (ntracks[ilayer]>95) break; //space for skipped clusters
1093 Bool_t changedet =kFALSE;
1094 if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1095 Int_t idetc=cl->GetDetectorIndex();
1097 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1098 // take into account misalignment (bring track to real detector plane)
1099 Double_t xTrOrig = currenttrack->GetX();
1100 currenttrack->PropagateTo(xTrOrig+cl->GetX(),0.,0.);
1101 // a first cut on track-cluster distance
1102 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1103 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1104 { // cluster not associated to track
1105 //printf("not ass\n");
1108 // bring track back to ideal detector plane
1109 currenttrack->PropagateTo(xTrOrig,0.,0.);
1110 } else { // have to move track to cluster's detector
1111 const AliITSdetector &detc=layer.GetDetector(idetc);
1112 // a first cut on track-cluster distance
1114 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1115 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1116 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1117 continue; // cluster not associated to track
1119 new (&backuptrack) AliITStrackMI(currenttrack2);
1121 currenttrack =¤ttrack2;
1122 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1123 new (currenttrack) AliITStrackMI(backuptrack);
1127 currenttrack->SetDetectorIndex(idetc);
1128 // Get again the budget to the primary vertex
1129 // for the current track being prolonged, if had to change detector
1130 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1133 // calculate track-clusters chi2
1134 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1136 //printf("chi2 %f max %f\n",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer));
1137 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1138 if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1139 if (ntracks[ilayer]>=100) continue;
1140 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1141 updatetrack->SetClIndex(ilayer,0);
1142 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1144 if (cl->GetQ()!=0) { // real cluster
1145 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1146 //printf("update failed\n");
1149 updatetrack->SetSampledEdx(cl->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
1150 modstatus = 1; // found
1151 } else { // virtual cluster in dead zone
1152 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1153 updatetrack->SetDeadZoneProbability(GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1154 modstatus = 7; // holes in z in SPD
1158 Float_t xlocnewdet,zlocnewdet;
1159 LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet); // local module coords
1160 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1162 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1164 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1166 // apply correction for material of the current layer
1167 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1169 if (constrain) { // apply vertex constrain
1170 updatetrack->SetConstrain(constrain);
1171 Bool_t isPrim = kTRUE;
1172 if (ilayer<4) { // check that it's close to the vertex
1173 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1174 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1175 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1176 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1177 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1179 if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1180 } //apply vertex constrain
1182 } // create new hypothesis
1183 //else printf("chi2 too large\n");
1184 } // loop over possible prolongations
1186 // allow one prolongation without clusters
1187 if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1188 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1189 // apply correction for material of the current layer
1190 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1191 vtrack->SetClIndex(ilayer,0);
1192 modstatus = 3; // skipped
1193 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1194 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1195 vtrack->IncrementNSkipped();
1199 // allow one prolongation without clusters for tracks with |tgl|>1.1
1200 if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) { //big theta - for low flux
1201 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1202 // apply correction for material of the current layer
1203 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1204 vtrack->SetClIndex(ilayer,0);
1205 modstatus = 3; // skipped
1206 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1207 vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1208 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1213 } // loop over tracks in layer ilayer+1
1215 //loop over track candidates for the current layer
1221 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1222 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1223 if (normalizedchi2[itrack] <
1224 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1228 if (constrain) { // constrain
1229 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1231 } else { // !constrain
1232 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1237 // sort tracks by increasing normalized chi2
1238 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1239 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1240 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1241 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1242 } // end loop over layers
1246 // Now select tracks to be kept
1248 Int_t max = constrain ? 20 : 5;
1250 // tracks that reach layer 0 (SPD inner)
1251 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1252 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1253 if (track.GetNumberOfClusters()<2) continue;
1254 if (!constrain && track.GetNormChi2(0) >
1255 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1258 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1261 // tracks that reach layer 1 (SPD outer)
1262 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1263 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1264 if (track.GetNumberOfClusters()<4) continue;
1265 if (!constrain && track.GetNormChi2(1) >
1266 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1267 if (constrain) track.IncrementNSkipped();
1269 track.SetD(0,track.GetD(GetX(),GetY()));
1270 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1271 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1272 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1275 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1278 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1280 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1281 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1282 if (track.GetNumberOfClusters()<3) continue;
1283 if (!constrain && track.GetNormChi2(2) >
1284 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1285 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1287 track.SetD(0,track.GetD(GetX(),GetY()));
1288 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1289 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1290 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1293 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1299 // register best track of each layer - important for V0 finder
1301 for (Int_t ilayer=0;ilayer<5;ilayer++){
1302 if (ntracks[ilayer]==0) continue;
1303 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1304 if (track.GetNumberOfClusters()<1) continue;
1305 CookLabel(&track,0);
1306 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1310 // update TPC V0 information
1312 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1313 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1314 for (Int_t i=0;i<3;i++){
1315 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1316 if (index==0) break;
1317 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1318 if (vertex->GetStatus()<0) continue; // rejected V0
1320 if (otrack->GetSign()>0) {
1321 vertex->SetIndex(0,esdindex);
1324 vertex->SetIndex(1,esdindex);
1326 //find nearest layer with track info
1327 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1328 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1329 Int_t nearest = nearestold;
1330 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1331 if (ntracks[nearest]==0){
1336 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1337 if (nearestold<5&&nearest<5){
1338 Bool_t accept = track.GetNormChi2(nearest)<10;
1340 if (track.GetSign()>0) {
1341 vertex->SetParamP(track);
1342 vertex->Update(fprimvertex);
1343 //vertex->SetIndex(0,track.fESDtrack->GetID());
1344 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1346 vertex->SetParamN(track);
1347 vertex->Update(fprimvertex);
1348 //vertex->SetIndex(1,track.fESDtrack->GetID());
1349 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1351 vertex->SetStatus(vertex->GetStatus()+1);
1353 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1360 //------------------------------------------------------------------------
1361 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1363 //--------------------------------------------------------------------
1366 return fgLayers[layer];
1368 //------------------------------------------------------------------------
1369 AliITStrackerMI::AliITSlayer::AliITSlayer():
1394 //--------------------------------------------------------------------
1395 //default AliITSlayer constructor
1396 //--------------------------------------------------------------------
1397 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1398 fClusterWeight[i]=0;
1399 fClusterTracks[0][i]=-1;
1400 fClusterTracks[1][i]=-1;
1401 fClusterTracks[2][i]=-1;
1402 fClusterTracks[3][i]=-1;
1405 //------------------------------------------------------------------------
1406 AliITStrackerMI::AliITSlayer::
1407 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1432 //--------------------------------------------------------------------
1433 //main AliITSlayer constructor
1434 //--------------------------------------------------------------------
1435 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1436 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1438 //------------------------------------------------------------------------
1439 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1441 fPhiOffset(layer.fPhiOffset),
1442 fNladders(layer.fNladders),
1443 fZOffset(layer.fZOffset),
1444 fNdetectors(layer.fNdetectors),
1445 fDetectors(layer.fDetectors),
1450 fClustersCs(layer.fClustersCs),
1451 fClusterIndexCs(layer.fClusterIndexCs),
1455 fCurrentSlice(layer.fCurrentSlice),
1462 fAccepted(layer.fAccepted),
1466 //------------------------------------------------------------------------
1467 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1468 //--------------------------------------------------------------------
1469 // AliITSlayer destructor
1470 //--------------------------------------------------------------------
1471 delete [] fDetectors;
1472 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1473 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1474 fClusterWeight[i]=0;
1475 fClusterTracks[0][i]=-1;
1476 fClusterTracks[1][i]=-1;
1477 fClusterTracks[2][i]=-1;
1478 fClusterTracks[3][i]=-1;
1481 //------------------------------------------------------------------------
1482 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1483 //--------------------------------------------------------------------
1484 // This function removes loaded clusters
1485 //--------------------------------------------------------------------
1486 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1487 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1488 fClusterWeight[i]=0;
1489 fClusterTracks[0][i]=-1;
1490 fClusterTracks[1][i]=-1;
1491 fClusterTracks[2][i]=-1;
1492 fClusterTracks[3][i]=-1;
1498 //------------------------------------------------------------------------
1499 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1500 //--------------------------------------------------------------------
1501 // This function reset weights of the clusters
1502 //--------------------------------------------------------------------
1503 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1504 fClusterWeight[i]=0;
1505 fClusterTracks[0][i]=-1;
1506 fClusterTracks[1][i]=-1;
1507 fClusterTracks[2][i]=-1;
1508 fClusterTracks[3][i]=-1;
1510 for (Int_t i=0; i<fN;i++) {
1511 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1512 if (cl&&cl->IsUsed()) cl->Use();
1516 //------------------------------------------------------------------------
1517 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1518 //--------------------------------------------------------------------
1519 // This function calculates the road defined by the cluster density
1520 //--------------------------------------------------------------------
1522 for (Int_t i=0; i<fN; i++) {
1523 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1525 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1527 //------------------------------------------------------------------------
1528 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1529 //--------------------------------------------------------------------
1530 //This function adds a cluster to this layer
1531 //--------------------------------------------------------------------
1532 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1533 ::Error("InsertCluster","Too many clusters !\n");
1539 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1540 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1541 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1542 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1543 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1547 //------------------------------------------------------------------------
1548 void AliITStrackerMI::AliITSlayer::SortClusters()
1553 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1554 Float_t *z = new Float_t[fN];
1555 Int_t * index = new Int_t[fN];
1557 for (Int_t i=0;i<fN;i++){
1558 z[i] = fClusters[i]->GetZ();
1560 TMath::Sort(fN,z,index,kFALSE);
1561 for (Int_t i=0;i<fN;i++){
1562 clusters[i] = fClusters[index[i]];
1565 for (Int_t i=0;i<fN;i++){
1566 fClusters[i] = clusters[i];
1567 fZ[i] = fClusters[i]->GetZ();
1568 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1569 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1570 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1580 for (Int_t i=0;i<fN;i++){
1581 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1582 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1583 fClusterIndex[i] = i;
1587 fDy5 = (fYB[1]-fYB[0])/5.;
1588 fDy10 = (fYB[1]-fYB[0])/10.;
1589 fDy20 = (fYB[1]-fYB[0])/20.;
1590 for (Int_t i=0;i<6;i++) fN5[i] =0;
1591 for (Int_t i=0;i<11;i++) fN10[i]=0;
1592 for (Int_t i=0;i<21;i++) fN20[i]=0;
1594 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;}
1595 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;}
1596 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;}
1599 for (Int_t i=0;i<fN;i++)
1600 for (Int_t irot=-1;irot<=1;irot++){
1601 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1603 for (Int_t slice=0; slice<6;slice++){
1604 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1605 fClusters5[slice][fN5[slice]] = fClusters[i];
1606 fY5[slice][fN5[slice]] = curY;
1607 fZ5[slice][fN5[slice]] = fZ[i];
1608 fClusterIndex5[slice][fN5[slice]]=i;
1613 for (Int_t slice=0; slice<11;slice++){
1614 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1615 fClusters10[slice][fN10[slice]] = fClusters[i];
1616 fY10[slice][fN10[slice]] = curY;
1617 fZ10[slice][fN10[slice]] = fZ[i];
1618 fClusterIndex10[slice][fN10[slice]]=i;
1623 for (Int_t slice=0; slice<21;slice++){
1624 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1625 fClusters20[slice][fN20[slice]] = fClusters[i];
1626 fY20[slice][fN20[slice]] = curY;
1627 fZ20[slice][fN20[slice]] = fZ[i];
1628 fClusterIndex20[slice][fN20[slice]]=i;
1635 // consistency check
1637 for (Int_t i=0;i<fN-1;i++){
1643 for (Int_t slice=0;slice<21;slice++)
1644 for (Int_t i=0;i<fN20[slice]-1;i++){
1645 if (fZ20[slice][i]>fZ20[slice][i+1]){
1652 //------------------------------------------------------------------------
1653 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1654 //--------------------------------------------------------------------
1655 // This function returns the index of the nearest cluster
1656 //--------------------------------------------------------------------
1659 if (fCurrentSlice<0) {
1668 if (ncl==0) return 0;
1669 Int_t b=0, e=ncl-1, m=(b+e)/2;
1670 for (; b<e; m=(b+e)/2) {
1671 // if (z > fClusters[m]->GetZ()) b=m+1;
1672 if (z > zcl[m]) b=m+1;
1677 //------------------------------------------------------------------------
1678 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 {
1679 //--------------------------------------------------------------------
1680 // This function computes the rectangular road for this track
1681 //--------------------------------------------------------------------
1684 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1685 // take into account the misalignment: propagate track to misaligned detector plane
1686 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1688 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1689 TMath::Sqrt(track->GetSigmaZ2() +
1690 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1691 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1692 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1693 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1694 TMath::Sqrt(track->GetSigmaY2() +
1695 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1696 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1697 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1699 // track at boundary between detectors, enlarge road
1700 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1701 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1702 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1703 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1704 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1705 Float_t tgl = TMath::Abs(track->GetTgl());
1706 if (tgl > 1.) tgl=1.;
1707 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1708 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1709 Float_t snp = TMath::Abs(track->GetSnp());
1710 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1711 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1714 // add to the road a term (up to 2-3 mm) to deal with misalignments
1715 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1716 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1718 Double_t r = fgLayers[ilayer].GetR();
1719 zmin = track->GetZ() - dz;
1720 zmax = track->GetZ() + dz;
1721 ymin = track->GetY() + r*det.GetPhi() - dy;
1722 ymax = track->GetY() + r*det.GetPhi() + dy;
1724 // bring track back to idead detector plane
1725 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1729 //------------------------------------------------------------------------
1730 void AliITStrackerMI::AliITSlayer::
1731 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1732 //--------------------------------------------------------------------
1733 // This function sets the "window"
1734 //--------------------------------------------------------------------
1736 Double_t circle=2*TMath::Pi()*fR;
1737 fYmin = ymin; fYmax =ymax;
1738 Float_t ymiddle = (fYmin+fYmax)*0.5;
1739 if (ymiddle<fYB[0]) {
1740 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1741 } else if (ymiddle>fYB[1]) {
1742 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1748 fClustersCs = fClusters;
1749 fClusterIndexCs = fClusterIndex;
1755 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1756 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1757 if (slice<0) slice=0;
1758 if (slice>20) slice=20;
1759 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1761 fCurrentSlice=slice;
1762 fClustersCs = fClusters20[fCurrentSlice];
1763 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1764 fYcs = fY20[fCurrentSlice];
1765 fZcs = fZ20[fCurrentSlice];
1766 fNcs = fN20[fCurrentSlice];
1771 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1772 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1773 if (slice<0) slice=0;
1774 if (slice>10) slice=10;
1775 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1777 fCurrentSlice=slice;
1778 fClustersCs = fClusters10[fCurrentSlice];
1779 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1780 fYcs = fY10[fCurrentSlice];
1781 fZcs = fZ10[fCurrentSlice];
1782 fNcs = fN10[fCurrentSlice];
1787 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1788 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1789 if (slice<0) slice=0;
1790 if (slice>5) slice=5;
1791 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1793 fCurrentSlice=slice;
1794 fClustersCs = fClusters5[fCurrentSlice];
1795 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1796 fYcs = fY5[fCurrentSlice];
1797 fZcs = fZ5[fCurrentSlice];
1798 fNcs = fN5[fCurrentSlice];
1802 fI=FindClusterIndex(zmin); fZmax=zmax;
1803 fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1809 //------------------------------------------------------------------------
1810 Int_t AliITStrackerMI::AliITSlayer::
1811 FindDetectorIndex(Double_t phi, Double_t z) const {
1812 //--------------------------------------------------------------------
1813 //This function finds the detector crossed by the track
1814 //--------------------------------------------------------------------
1816 if (fZOffset<0) // old geometry
1817 dphi = -(phi-fPhiOffset);
1818 else // new geometry
1819 dphi = phi-fPhiOffset;
1822 if (dphi < 0) dphi += 2*TMath::Pi();
1823 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1824 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1825 if (np>=fNladders) np-=fNladders;
1826 if (np<0) np+=fNladders;
1829 Double_t dz=fZOffset-z;
1830 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1831 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1832 if (nz>=fNdetectors) return -1;
1833 if (nz<0) return -1;
1835 // ad hoc correction for 3rd ladder of SDD inner layer,
1836 // which is reversed (rotated by pi around local y)
1837 // this correction is OK only from AliITSv11Hybrid onwards
1838 if (GetR()>12. && GetR()<20.) { // SDD inner
1839 if(np==2) { // 3rd ladder
1840 nz = (fNdetectors-1) - nz;
1843 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1846 return np*fNdetectors + nz;
1848 //------------------------------------------------------------------------
1849 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1851 //--------------------------------------------------------------------
1852 // This function returns clusters within the "window"
1853 //--------------------------------------------------------------------
1855 if (fCurrentSlice<0) {
1856 Double_t rpi2 = 2.*fR*TMath::Pi();
1857 for (Int_t i=fI; i<fImax; i++) {
1859 if (fYmax<y) y -= rpi2;
1860 if (fYmin>y) y += rpi2;
1861 if (y<fYmin) continue;
1862 if (y>fYmax) continue;
1863 if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1866 return fClusters[i];
1869 for (Int_t i=fI; i<fImax; i++) {
1870 if (fYcs[i]<fYmin) continue;
1871 if (fYcs[i]>fYmax) continue;
1872 if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1873 ci=fClusterIndexCs[i];
1875 return fClustersCs[i];
1880 //------------------------------------------------------------------------
1881 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1883 //--------------------------------------------------------------------
1884 // This function returns the layer thickness at this point (units X0)
1885 //--------------------------------------------------------------------
1887 x0=AliITSRecoParam::GetX0Air();
1888 if (43<fR&&fR<45) { //SSD2
1891 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1892 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1893 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1894 for (Int_t i=0; i<12; i++) {
1895 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1896 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1900 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1901 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1905 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1906 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1909 if (37<fR&&fR<41) { //SSD1
1912 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1913 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1914 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1915 for (Int_t i=0; i<11; i++) {
1916 if (TMath::Abs(z-3.9*i)<0.15) {
1917 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1921 if (TMath::Abs(z+3.9*i)<0.15) {
1922 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1926 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1927 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1930 if (13<fR&&fR<26) { //SDD
1933 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1935 if (TMath::Abs(y-1.80)<0.55) {
1937 for (Int_t j=0; j<20; j++) {
1938 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1939 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1942 if (TMath::Abs(y+1.80)<0.55) {
1944 for (Int_t j=0; j<20; j++) {
1945 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1946 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1950 for (Int_t i=0; i<4; i++) {
1951 if (TMath::Abs(z-7.3*i)<0.60) {
1953 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1956 if (TMath::Abs(z+7.3*i)<0.60) {
1958 if (TMath::Abs(y-0.00)>3.30) d+=dd;
1963 if (6<fR&&fR<8) { //SPD2
1964 Double_t dd=0.0063; x0=21.5;
1966 if (TMath::Abs(y-3.08)>0.5) d+=dd;
1967 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
1969 if (3<fR&&fR<5) { //SPD1
1970 Double_t dd=0.0063; x0=21.5;
1972 if (TMath::Abs(y+0.21)>0.6) d+=dd;
1973 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
1978 //------------------------------------------------------------------------
1979 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
1981 fRmisal(det.fRmisal),
1983 fSinPhi(det.fSinPhi),
1984 fCosPhi(det.fCosPhi),
1990 fNChips(det.fNChips),
1991 fChipIsBad(det.fChipIsBad)
1995 //------------------------------------------------------------------------
1996 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
1997 AliITSDetTypeRec *detTypeRec)
1999 //--------------------------------------------------------------------
2000 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2001 //--------------------------------------------------------------------
2003 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2004 // while in the tracker they start from 0 for each layer
2005 for(Int_t il=0; il<ilayer; il++)
2006 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2009 if (ilayer==0 || ilayer==1) { // ---------- SPD
2011 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2013 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2016 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2020 // Get calibration from AliITSDetTypeRec
2021 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2022 AliITSCalibration *calibSPDdead = 0;
2023 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2024 if (calib->IsBad() ||
2025 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2028 printf("lay %d bad %d\n",ilayer,idet);
2031 // Get segmentation from AliITSDetTypeRec
2032 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2034 // Read info about bad chips
2035 fNChips = segm->GetMaximumChipIndex()+1;
2036 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2037 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2038 fChipIsBad = new Bool_t[fNChips];
2039 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2040 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2041 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2046 //------------------------------------------------------------------------
2047 Double_t AliITStrackerMI::GetEffectiveThickness()
2049 //--------------------------------------------------------------------
2050 // Returns the thickness between the current layer and the vertex (units X0)
2051 //--------------------------------------------------------------------
2054 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2055 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2056 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2060 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2061 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2065 Double_t xn=fgLayers[fI].GetR();
2066 for (Int_t i=0; i<fI; i++) {
2067 Double_t xi=fgLayers[i].GetR();
2068 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2074 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2075 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2078 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2079 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2083 //------------------------------------------------------------------------
2084 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2085 //-------------------------------------------------------------------
2086 // This function returns number of clusters within the "window"
2087 //--------------------------------------------------------------------
2089 for (Int_t i=fI; i<fN; i++) {
2090 const AliITSRecPoint *c=fClusters[i];
2091 if (c->GetZ() > fZmax) break;
2092 if (c->IsUsed()) continue;
2093 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2094 Double_t y=fR*det.GetPhi() + c->GetY();
2096 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2097 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2099 if (y<fYmin) continue;
2100 if (y>fYmax) continue;
2105 //------------------------------------------------------------------------
2106 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2107 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2109 //--------------------------------------------------------------------
2110 // This function refits the track "track" at the position "x" using
2111 // the clusters from "clusters"
2112 // If "extra"==kTRUE,
2113 // the clusters from overlapped modules get attached to "track"
2114 // If "planeff"==kTRUE,
2115 // special approach for plane efficiency evaluation is applyed
2116 //--------------------------------------------------------------------
2118 Int_t index[AliITSgeomTGeo::kNLayers];
2120 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2121 Int_t nc=clusters->GetNumberOfClusters();
2122 for (k=0; k<nc; k++) {
2123 Int_t idx=clusters->GetClusterIndex(k);
2124 Int_t ilayer=(idx&0xf0000000)>>28;
2128 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2130 //------------------------------------------------------------------------
2131 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2132 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2134 //--------------------------------------------------------------------
2135 // This function refits the track "track" at the position "x" using
2136 // the clusters from array
2137 // If "extra"==kTRUE,
2138 // the clusters from overlapped modules get attached to "track"
2139 // If "planeff"==kTRUE,
2140 // special approach for plane efficiency evaluation is applyed
2141 //--------------------------------------------------------------------
2142 Int_t index[AliITSgeomTGeo::kNLayers];
2144 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2146 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2147 index[k]=clusters[k];
2150 // special for cosmics: check which the innermost layer crossed
2152 Int_t innermostlayer=5;
2153 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2154 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2155 if(drphi < fgLayers[innermostlayer].GetR()) break;
2157 //printf(" drphi %f innermost %d\n",drphi,innermostlayer);
2159 Int_t modstatus=1; // found
2161 Int_t from, to, step;
2162 if (xx > track->GetX()) {
2163 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2166 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2169 TString dir = (step>0 ? "outward" : "inward");
2171 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2172 AliITSlayer &layer=fgLayers[ilayer];
2173 Double_t r=layer.GetR();
2174 if (step<0 && xx>r) break;
2176 // material between SSD and SDD, SDD and SPD
2177 Double_t hI=ilayer-0.5*step;
2178 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2179 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2180 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2181 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2183 // remember old position [SR, GSI 18.02.2003]
2184 Double_t oldX=0., oldY=0., oldZ=0.;
2185 if (track->IsStartedTimeIntegral() && step==1) {
2186 track->GetGlobalXYZat(track->GetX(),oldX,oldY,oldZ);
2190 Double_t oldGlobXYZ[3];
2191 track->GetXYZ(oldGlobXYZ);
2194 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2196 Int_t idet=layer.FindDetectorIndex(phi,z);
2198 // check if we allow a prolongation without point for large-eta tracks
2199 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2201 // propagate to the layer radius
2202 Double_t xToGo; track->GetLocalXat(r,xToGo);
2203 track->AliExternalTrackParam::PropagateTo(xToGo,GetBz());
2204 // apply correction for material of the current layer
2205 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2206 modstatus = 4; // out in z
2207 LocalModuleCoord(ilayer,idet,track,xloc,zloc); // local module coords
2208 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2209 // track time update [SR, GSI 17.02.2003]
2210 if (track->IsStartedTimeIntegral() && step==1) {
2211 Double_t newX, newY, newZ;
2212 track->GetGlobalXYZat(track->GetX(),newX,newY,newZ);
2213 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2214 (oldZ-newZ)*(oldZ-newZ);
2215 track->AddTimeStep(TMath::Sqrt(dL2));
2220 if (idet<0) return kFALSE;
2222 const AliITSdetector &det=layer.GetDetector(idet);
2223 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2225 track->SetDetectorIndex(idet);
2226 LocalModuleCoord(ilayer,idet,track,xloc,zloc); // local module coords
2228 Double_t dz,zmin,zmax,dy,ymin,ymax;
2230 const AliITSRecPoint *clAcc=0;
2231 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2233 Int_t idx=index[ilayer];
2234 if (idx>=0) { // cluster in this layer
2235 modstatus = 6; // no refit
2236 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2238 if (idet != cl->GetDetectorIndex()) {
2239 idet=cl->GetDetectorIndex();
2240 const AliITSdetector &detc=layer.GetDetector(idet);
2241 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2242 track->SetDetectorIndex(idet);
2243 LocalModuleCoord(ilayer,idet,track,xloc,zloc); // local module coords
2245 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2246 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2250 modstatus = 1; // found
2255 } else { // no cluster in this layer
2257 modstatus = 3; // skipped
2258 // Plane Eff determination:
2259 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2260 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2261 UseTrackForPlaneEff(track,ilayer);
2264 modstatus = 5; // no cls in road
2266 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2267 dz = 0.5*(zmax-zmin);
2268 dy = 0.5*(ymax-ymin);
2269 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2270 if (dead==1) modstatus = 7; // holes in z in SPD
2271 if (dead==2 || dead==3) modstatus = 2; // dead from OCDB
2276 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2277 track->SetSampledEdx(clAcc->GetQ(),track->GetNumberOfClusters()-1);
2279 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2282 if (extra) { // search for extra clusters in overlapped modules
2283 AliITStrackV2 tmp(*track);
2284 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2285 layer.SelectClusters(zmin,zmax,ymin,ymax);
2287 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2289 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2290 Double_t tolerance=0.1;
2291 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2292 // only clusters in another module! (overlaps)
2293 idetExtra = clExtra->GetDetectorIndex();
2294 if (idet == idetExtra) continue;
2296 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2298 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2299 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2300 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2301 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2303 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2304 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2307 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2308 track->SetExtraModule(ilayer,idetExtra);
2310 } // end search for extra clusters in overlapped modules
2312 // Correct for material of the current layer
2313 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2315 // track time update [SR, GSI 17.02.2003]
2316 if (track->IsStartedTimeIntegral() && step==1) {
2317 Double_t newX, newY, newZ;
2318 track->GetGlobalXYZat(track->GetX(),newX,newY,newZ);
2319 Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
2320 (oldZ-newZ)*(oldZ-newZ);
2321 track->AddTimeStep(TMath::Sqrt(dL2));
2325 } // end loop on layers
2327 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2331 //------------------------------------------------------------------------
2332 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2335 // calculate normalized chi2
2336 // return NormalizedChi2(track,0);
2339 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2340 // track->fdEdxMismatch=0;
2341 Float_t dedxmismatch =0;
2342 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2344 for (Int_t i = 0;i<6;i++){
2345 if (track->GetClIndex(i)>0){
2346 Float_t cerry, cerrz;
2347 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2349 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2352 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2353 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2354 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2356 cchi2+=(0.5-ratio)*10.;
2357 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2358 dedxmismatch+=(0.5-ratio)*10.;
2362 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2363 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2364 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2365 if (i<2) chi2+=2*cl->GetDeltaProbability();
2371 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2372 track->SetdEdxMismatch(dedxmismatch);
2376 for (Int_t i = 0;i<4;i++){
2377 if (track->GetClIndex(i)>0){
2378 Float_t cerry, cerrz;
2379 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2380 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2383 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2384 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2388 for (Int_t i = 4;i<6;i++){
2389 if (track->GetClIndex(i)>0){
2390 Float_t cerry, cerrz;
2391 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2392 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2395 Float_t cerryb, cerrzb;
2396 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2397 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2400 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2401 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2406 if (track->GetESDtrack()->GetTPCsignal()>85){
2407 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2409 chi2+=(0.5-ratio)*5.;
2412 chi2+=(ratio-2.0)*3;
2416 Double_t match = TMath::Sqrt(track->GetChi22());
2417 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2418 if (!track->GetConstrain()) {
2419 if (track->GetNumberOfClusters()>2) {
2420 match/=track->GetNumberOfClusters()-2.;
2425 if (match<0) match=0;
2426 Float_t deadzonefactor = (track->GetNDeadZone()>0) ? 3*(1.1-track->GetDeadZoneProbability()):0.;
2427 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2428 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2429 1./(1.+track->GetNSkipped()));
2433 //------------------------------------------------------------------------
2434 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
2437 // return matching chi2 between two tracks
2438 AliITStrackMI track3(*track2);
2439 track3.Propagate(track1->GetAlpha(),track1->GetX());
2441 vec(0,0)=track1->GetY() - track3.GetY();
2442 vec(1,0)=track1->GetZ() - track3.GetZ();
2443 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2444 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2445 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2448 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2449 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2450 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2451 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2452 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2454 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2455 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2456 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2457 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2459 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2460 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2461 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2463 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2464 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2466 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2469 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2470 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2473 //------------------------------------------------------------------------
2474 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr)
2477 // return probability that given point (characterized by z position and error)
2478 // is in SPD dead zone
2480 Double_t probability = 0.;
2481 Double_t absz = TMath::Abs(zpos);
2482 Double_t nearestz = (absz<2.) ? 0.5*(fSPDdetzcentre[1]+fSPDdetzcentre[2]) :
2483 0.5*(fSPDdetzcentre[2]+fSPDdetzcentre[3]);
2484 if (TMath::Abs(absz-nearestz)>0.25+3.*zerr) return probability;
2485 Double_t zmin, zmax;
2486 if (zpos<-6.) { // dead zone at z = -7
2487 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2488 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2489 } else if (zpos>6.) { // dead zone at z = +7
2490 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2491 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2492 } else if (absz<2.) { // dead zone at z = 0
2493 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2494 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2499 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2501 probability = 0.5*( TMath::Erf((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2502 TMath::Erf((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2505 //------------------------------------------------------------------------
2506 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
2509 // calculate normalized chi2
2511 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2513 for (Int_t i = 0;i<6;i++){
2514 if (TMath::Abs(track->GetDy(i))>0){
2515 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2516 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2519 else{chi2[i]=10000;}
2522 TMath::Sort(6,chi2,index,kFALSE);
2523 Float_t max = float(ncl)*fac-1.;
2524 Float_t sumchi=0, sumweight=0;
2525 for (Int_t i=0;i<max+1;i++){
2526 Float_t weight = (i<max)?1.:(max+1.-i);
2527 sumchi+=weight*chi2[index[i]];
2530 Double_t normchi2 = sumchi/sumweight;
2533 //------------------------------------------------------------------------
2534 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
2537 // calculate normalized chi2
2538 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2541 for (Int_t i=0;i<6;i++){
2542 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2543 Double_t sy1 = forwardtrack->GetSigmaY(i);
2544 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2545 Double_t sy2 = backtrack->GetSigmaY(i);
2546 Double_t sz2 = backtrack->GetSigmaZ(i);
2547 if (i<2){ sy2=1000.;sz2=1000;}
2549 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2550 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2552 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2553 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2555 res+= nz0*nz0+ny0*ny0;
2558 if (npoints>1) return
2559 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2560 //2*forwardtrack->fNUsed+
2561 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2562 1./(1.+forwardtrack->GetNSkipped()));
2565 //------------------------------------------------------------------------
2566 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2567 //--------------------------------------------------------------------
2568 // Return pointer to a given cluster
2569 //--------------------------------------------------------------------
2570 Int_t l=(index & 0xf0000000) >> 28;
2571 Int_t c=(index & 0x0fffffff) >> 00;
2572 return fgLayers[l].GetWeight(c);
2574 //------------------------------------------------------------------------
2575 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
2577 //---------------------------------------------
2578 // register track to the list
2580 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2583 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2584 Int_t index = track->GetClusterIndex(icluster);
2585 Int_t l=(index & 0xf0000000) >> 28;
2586 Int_t c=(index & 0x0fffffff) >> 00;
2587 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2588 for (Int_t itrack=0;itrack<4;itrack++){
2589 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2590 fgLayers[l].SetClusterTracks(itrack,c,id);
2596 //------------------------------------------------------------------------
2597 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
2599 //---------------------------------------------
2600 // unregister track from the list
2601 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2602 Int_t index = track->GetClusterIndex(icluster);
2603 Int_t l=(index & 0xf0000000) >> 28;
2604 Int_t c=(index & 0x0fffffff) >> 00;
2605 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2606 for (Int_t itrack=0;itrack<4;itrack++){
2607 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2608 fgLayers[l].SetClusterTracks(itrack,c,-1);
2613 //------------------------------------------------------------------------
2614 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2616 //-------------------------------------------------------------
2617 //get number of shared clusters
2618 //-------------------------------------------------------------
2620 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2621 // mean number of clusters
2622 Float_t *ny = GetNy(id), *nz = GetNz(id);
2625 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2626 Int_t index = track->GetClusterIndex(icluster);
2627 Int_t l=(index & 0xf0000000) >> 28;
2628 Int_t c=(index & 0x0fffffff) >> 00;
2629 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2631 printf("problem\n");
2633 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2637 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2638 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2639 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2641 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2644 deltan = (cl->GetNz()-nz[l]);
2646 if (deltan>2.0) continue; // extended - highly probable shared cluster
2647 weight = 2./TMath::Max(3.+deltan,2.);
2649 for (Int_t itrack=0;itrack<4;itrack++){
2650 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2652 clist[l] = (AliITSRecPoint*)GetCluster(index);
2658 track->SetNUsed(shared);
2661 //------------------------------------------------------------------------
2662 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2665 // find first shared track
2667 // mean number of clusters
2668 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2670 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2671 Int_t sharedtrack=100000;
2672 Int_t tracks[24],trackindex=0;
2673 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2675 for (Int_t icluster=0;icluster<6;icluster++){
2676 if (clusterlist[icluster]<0) continue;
2677 Int_t index = clusterlist[icluster];
2678 Int_t l=(index & 0xf0000000) >> 28;
2679 Int_t c=(index & 0x0fffffff) >> 00;
2681 printf("problem\n");
2683 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2684 //if (l>3) continue;
2685 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2688 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2689 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2690 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2692 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2695 deltan = (cl->GetNz()-nz[l]);
2697 if (deltan>2.0) continue; // extended - highly probable shared cluster
2699 for (Int_t itrack=3;itrack>=0;itrack--){
2700 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2701 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2702 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2707 if (trackindex==0) return -1;
2709 sharedtrack = tracks[0];
2711 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2714 Int_t tracks2[24], cluster[24];
2715 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2718 for (Int_t i=0;i<trackindex;i++){
2719 if (tracks[i]<0) continue;
2720 tracks2[index] = tracks[i];
2722 for (Int_t j=i+1;j<trackindex;j++){
2723 if (tracks[j]<0) continue;
2724 if (tracks[j]==tracks[i]){
2732 for (Int_t i=0;i<index;i++){
2733 if (cluster[index]>max) {
2734 sharedtrack=tracks2[index];
2741 if (sharedtrack>=100000) return -1;
2743 // make list of overlaps
2745 for (Int_t icluster=0;icluster<6;icluster++){
2746 if (clusterlist[icluster]<0) continue;
2747 Int_t index = clusterlist[icluster];
2748 Int_t l=(index & 0xf0000000) >> 28;
2749 Int_t c=(index & 0x0fffffff) >> 00;
2750 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2751 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2753 if (cl->GetNy()>2) continue;
2754 if (cl->GetNz()>2) continue;
2757 if (cl->GetNy()>3) continue;
2758 if (cl->GetNz()>3) continue;
2761 for (Int_t itrack=3;itrack>=0;itrack--){
2762 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2763 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2771 //------------------------------------------------------------------------
2772 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2774 // try to find track hypothesys without conflicts
2775 // with minimal chi2;
2776 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2777 Int_t entries1 = arr1->GetEntriesFast();
2778 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2779 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2780 Int_t entries2 = arr2->GetEntriesFast();
2781 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2783 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2784 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2785 if (track10->Pt()>0.5+track20->Pt()) return track10;
2787 for (Int_t itrack=0;itrack<entries1;itrack++){
2788 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2789 UnRegisterClusterTracks(track,trackID1);
2792 for (Int_t itrack=0;itrack<entries2;itrack++){
2793 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2794 UnRegisterClusterTracks(track,trackID2);
2798 Float_t maxconflicts=6;
2799 Double_t maxchi2 =1000.;
2801 // get weight of hypothesys - using best hypothesy
2804 Int_t list1[6],list2[6];
2805 AliITSRecPoint *clist1[6], *clist2[6] ;
2806 RegisterClusterTracks(track10,trackID1);
2807 RegisterClusterTracks(track20,trackID2);
2808 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2809 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2810 UnRegisterClusterTracks(track10,trackID1);
2811 UnRegisterClusterTracks(track20,trackID2);
2814 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2815 Float_t nerry[6],nerrz[6];
2816 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2817 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2818 for (Int_t i=0;i<6;i++){
2819 if ( (erry1[i]>0) && (erry2[i]>0)) {
2820 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2821 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2823 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2824 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2826 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2827 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2828 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2831 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2832 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2833 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2841 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2842 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2843 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2844 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
2846 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
2847 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
2848 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
2850 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
2851 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
2852 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
2855 Double_t sumw = w1+w2;
2859 w1 = (d2+0.5)/(d1+d2+1.);
2860 w2 = (d1+0.5)/(d1+d2+1.);
2862 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
2863 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
2865 // get pair of "best" hypothesys
2867 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
2868 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
2870 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2871 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2872 //if (track1->fFakeRatio>0) continue;
2873 RegisterClusterTracks(track1,trackID1);
2874 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2875 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2877 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2878 //if (track2->fFakeRatio>0) continue;
2880 RegisterClusterTracks(track2,trackID2);
2881 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2882 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2883 UnRegisterClusterTracks(track2,trackID2);
2885 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
2886 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
2887 if (nskipped>0.5) continue;
2889 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2890 if (conflict1+1<cconflict1) continue;
2891 if (conflict2+1<cconflict2) continue;
2895 for (Int_t i=0;i<6;i++){
2897 Float_t c1 =0.; // conflict coeficients
2899 if (clist1[i]&&clist2[i]){
2902 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2905 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2907 c1 = 2./TMath::Max(3.+deltan,2.);
2908 c2 = 2./TMath::Max(3.+deltan,2.);
2914 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2917 deltan = (clist1[i]->GetNz()-nz1[i]);
2919 c1 = 2./TMath::Max(3.+deltan,2.);
2926 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2929 deltan = (clist2[i]->GetNz()-nz2[i]);
2931 c2 = 2./TMath::Max(3.+deltan,2.);
2937 if (TMath::Abs(track1->GetDy(i))>0.) {
2938 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
2939 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
2940 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2941 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
2943 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
2946 if (TMath::Abs(track2->GetDy(i))>0.) {
2947 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
2948 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
2949 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2950 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2953 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
2955 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2956 if (chi21>0) sum+=w1;
2957 if (chi22>0) sum+=w2;
2960 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
2961 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
2962 Double_t normchi2 = 2*conflict+sumchi2/sum;
2963 if ( normchi2 <maxchi2 ){
2966 maxconflicts = conflict;
2970 UnRegisterClusterTracks(track1,trackID1);
2973 // if (maxconflicts<4 && maxchi2<th0){
2974 if (maxchi2<th0*2.){
2975 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
2976 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2977 track1->SetChi2MIP(5,maxconflicts);
2978 track1->SetChi2MIP(6,maxchi2);
2979 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
2980 // track1->UpdateESDtrack(AliESDtrack::kITSin);
2981 track1->SetChi2MIP(8,index1);
2982 fBestTrackIndex[trackID1] =index1;
2983 UpdateESDtrack(track1, AliESDtrack::kITSin);
2985 else if (track10->GetChi2MIP(0)<th1){
2986 track10->SetChi2MIP(5,maxconflicts);
2987 track10->SetChi2MIP(6,maxchi2);
2988 // track10->UpdateESDtrack(AliESDtrack::kITSin);
2989 UpdateESDtrack(track10,AliESDtrack::kITSin);
2992 for (Int_t itrack=0;itrack<entries1;itrack++){
2993 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2994 UnRegisterClusterTracks(track,trackID1);
2997 for (Int_t itrack=0;itrack<entries2;itrack++){
2998 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2999 UnRegisterClusterTracks(track,trackID2);
3002 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3003 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3004 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3005 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3006 RegisterClusterTracks(track10,trackID1);
3008 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3009 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3010 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3011 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3012 RegisterClusterTracks(track20,trackID2);
3017 //------------------------------------------------------------------------
3018 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3019 //--------------------------------------------------------------------
3020 // This function marks clusters assigned to the track
3021 //--------------------------------------------------------------------
3022 AliTracker::UseClusters(t,from);
3024 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3025 //if (c->GetQ()>2) c->Use();
3026 if (c->GetSigmaZ2()>0.1) c->Use();
3027 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3028 //if (c->GetQ()>2) c->Use();
3029 if (c->GetSigmaZ2()>0.1) c->Use();
3032 //------------------------------------------------------------------------
3033 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3035 //------------------------------------------------------------------
3036 // add track to the list of hypothesys
3037 //------------------------------------------------------------------
3039 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3040 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3042 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3044 array = new TObjArray(10);
3045 fTrackHypothesys.AddAt(array,esdindex);
3047 array->AddLast(track);
3049 //------------------------------------------------------------------------
3050 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3052 //-------------------------------------------------------------------
3053 // compress array of track hypothesys
3054 // keep only maxsize best hypothesys
3055 //-------------------------------------------------------------------
3056 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3057 if (! (fTrackHypothesys.At(esdindex)) ) return;
3058 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3059 Int_t entries = array->GetEntriesFast();
3061 //- find preliminary besttrack as a reference
3062 Float_t minchi2=10000;
3064 AliITStrackMI * besttrack=0;
3065 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3066 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3067 if (!track) continue;
3068 Float_t chi2 = NormalizedChi2(track,0);
3070 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3071 track->SetLabel(tpcLabel);
3073 track->SetFakeRatio(1.);
3074 CookLabel(track,0.); //For comparison only
3076 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3077 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3078 if (track->GetNumberOfClusters()<maxn) continue;
3079 maxn = track->GetNumberOfClusters();
3086 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3087 delete array->RemoveAt(itrack);
3091 if (!besttrack) return;
3094 //take errors of best track as a reference
3095 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3096 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3097 for (Int_t j=0;j<6;j++) {
3098 if (besttrack->GetClIndex(j)>0){
3099 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3100 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3101 ny[j] = besttrack->GetNy(j);
3102 nz[j] = besttrack->GetNz(j);
3106 // calculate normalized chi2
3108 Float_t * chi2 = new Float_t[entries];
3109 Int_t * index = new Int_t[entries];
3110 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3111 for (Int_t itrack=0;itrack<entries;itrack++){
3112 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3114 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3115 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3116 chi2[itrack] = track->GetChi2MIP(0);
3118 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3119 delete array->RemoveAt(itrack);
3125 TMath::Sort(entries,chi2,index,kFALSE);
3126 besttrack = (AliITStrackMI*)array->At(index[0]);
3127 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3128 for (Int_t j=0;j<6;j++){
3129 if (besttrack->GetClIndex(j)>0){
3130 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3131 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3132 ny[j] = besttrack->GetNy(j);
3133 nz[j] = besttrack->GetNz(j);
3138 // calculate one more time with updated normalized errors
3139 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3140 for (Int_t itrack=0;itrack<entries;itrack++){
3141 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3143 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3144 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3145 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3148 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3149 delete array->RemoveAt(itrack);
3154 entries = array->GetEntriesFast();
3158 TObjArray * newarray = new TObjArray();
3159 TMath::Sort(entries,chi2,index,kFALSE);
3160 besttrack = (AliITStrackMI*)array->At(index[0]);
3163 for (Int_t j=0;j<6;j++){
3164 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3165 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3166 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3167 ny[j] = besttrack->GetNy(j);
3168 nz[j] = besttrack->GetNz(j);
3171 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3172 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3173 Float_t minn = besttrack->GetNumberOfClusters()-3;
3175 for (Int_t i=0;i<entries;i++){
3176 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3177 if (!track) continue;
3178 if (accepted>maxcut) break;
3179 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3180 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3181 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3182 delete array->RemoveAt(index[i]);
3186 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3187 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3188 if (!shortbest) accepted++;
3190 newarray->AddLast(array->RemoveAt(index[i]));
3191 for (Int_t j=0;j<6;j++){
3193 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3194 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3195 ny[j] = track->GetNy(j);
3196 nz[j] = track->GetNz(j);
3201 delete array->RemoveAt(index[i]);
3205 delete fTrackHypothesys.RemoveAt(esdindex);
3206 fTrackHypothesys.AddAt(newarray,esdindex);
3210 delete fTrackHypothesys.RemoveAt(esdindex);
3216 //------------------------------------------------------------------------
3217 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3219 //-------------------------------------------------------------
3220 // try to find best hypothesy
3221 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3222 //-------------------------------------------------------------
3223 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3224 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3225 if (!array) return 0;
3226 Int_t entries = array->GetEntriesFast();
3227 if (!entries) return 0;
3228 Float_t minchi2 = 100000;
3229 AliITStrackMI * besttrack=0;
3231 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3232 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3233 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3234 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3236 for (Int_t i=0;i<entries;i++){
3237 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3238 if (!track) continue;
3239 Float_t sigmarfi,sigmaz;
3240 GetDCASigma(track,sigmarfi,sigmaz);
3241 track->SetDnorm(0,sigmarfi);
3242 track->SetDnorm(1,sigmaz);
3244 track->SetChi2MIP(1,1000000);
3245 track->SetChi2MIP(2,1000000);
3246 track->SetChi2MIP(3,1000000);
3249 backtrack = new(backtrack) AliITStrackMI(*track);
3250 if (track->GetConstrain()) {
3251 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3252 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3253 backtrack->ResetCovariance(10.);
3255 backtrack->ResetCovariance(10.);
3257 backtrack->ResetClusters();
3259 Double_t x = original->GetX();
3260 if (!RefitAt(x,backtrack,track)) continue;
3262 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3263 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3264 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3265 track->SetChi22(GetMatchingChi2(backtrack,original));
3267 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3268 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3269 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3272 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3274 //forward track - without constraint
3275 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3276 forwardtrack->ResetClusters();
3278 RefitAt(x,forwardtrack,track);
3279 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3280 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3281 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3283 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3284 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3285 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3286 forwardtrack->SetD(0,track->GetD(0));
3287 forwardtrack->SetD(1,track->GetD(1));
3290 AliITSRecPoint* clist[6];
3291 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3292 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3295 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3296 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3297 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3298 track->SetChi2MIP(3,1000);
3301 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3303 for (Int_t ichi=0;ichi<5;ichi++){
3304 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3306 if (chi2 < minchi2){
3307 //besttrack = new AliITStrackMI(*forwardtrack);
3309 besttrack->SetLabel(track->GetLabel());
3310 besttrack->SetFakeRatio(track->GetFakeRatio());
3312 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3313 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3314 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3318 delete forwardtrack;
3320 for (Int_t i=0;i<entries;i++){
3321 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3323 if (!track) continue;
3325 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3326 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3327 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3328 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3329 delete array->RemoveAt(i);
3339 SortTrackHypothesys(esdindex,checkmax,1);
3341 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3342 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3343 besttrack = (AliITStrackMI*)array->At(0);
3344 if (!besttrack) return 0;
3345 besttrack->SetChi2MIP(8,0);
3346 fBestTrackIndex[esdindex]=0;
3347 entries = array->GetEntriesFast();
3348 AliITStrackMI *longtrack =0;
3350 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3351 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3352 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3353 if (!track->GetConstrain()) continue;
3354 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3355 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3356 if (track->GetChi2MIP(0)>4.) continue;
3357 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3360 //if (longtrack) besttrack=longtrack;
3363 AliITSRecPoint * clist[6];
3364 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3365 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3366 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3367 RegisterClusterTracks(besttrack,esdindex);
3374 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3375 if (sharedtrack>=0){
3377 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3379 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3385 if (shared>2.5) return 0;
3386 if (shared>1.0) return besttrack;
3388 // Don't sign clusters if not gold track
3390 if (!besttrack->IsGoldPrimary()) return besttrack;
3391 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3393 if (fConstraint[fPass]){
3397 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3398 for (Int_t i=0;i<6;i++){
3399 Int_t index = besttrack->GetClIndex(i);
3400 if (index<=0) continue;
3401 Int_t ilayer = (index & 0xf0000000) >> 28;
3402 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3403 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3405 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3406 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3407 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3408 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3409 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3410 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3412 Bool_t cansign = kTRUE;
3413 for (Int_t itrack=0;itrack<entries; itrack++){
3414 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3415 if (!track) continue;
3416 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3417 if ( (track->GetClIndex(ilayer)>0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3423 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3424 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3425 if (!c->IsUsed()) c->Use();
3431 //------------------------------------------------------------------------
3432 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3435 // get "best" hypothesys
3438 Int_t nentries = itsTracks.GetEntriesFast();
3439 for (Int_t i=0;i<nentries;i++){
3440 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3441 if (!track) continue;
3442 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3443 if (!array) continue;
3444 if (array->GetEntriesFast()<=0) continue;
3446 AliITStrackMI* longtrack=0;
3448 Float_t maxchi2=1000;
3449 for (Int_t j=0;j<array->GetEntriesFast();j++){
3450 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3451 if (!trackHyp) continue;
3452 if (trackHyp->GetGoldV0()) {
3453 longtrack = trackHyp; //gold V0 track taken
3456 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3457 Float_t chi2 = trackHyp->GetChi2MIP(0);
3459 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3461 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3463 if (chi2 > maxchi2) continue;
3464 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3471 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3472 if (!longtrack) {longtrack = besttrack;}
3473 else besttrack= longtrack;
3477 AliITSRecPoint * clist[6];
3478 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3480 track->SetNUsed(shared);
3481 track->SetNSkipped(besttrack->GetNSkipped());
3482 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3484 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3488 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3489 //if (sharedtrack==-1) sharedtrack=0;
3490 if (sharedtrack>=0) {
3491 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3494 if (besttrack&&fAfterV0) {
3495 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3497 if (besttrack&&fConstraint[fPass])
3498 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3499 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3500 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3501 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3507 //------------------------------------------------------------------------
3508 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3509 //--------------------------------------------------------------------
3510 //This function "cooks" a track label. If label<0, this track is fake.
3511 //--------------------------------------------------------------------
3514 if ( track->GetESDtrack()) tpcLabel = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3516 track->SetChi2MIP(9,0);
3518 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3519 Int_t cindex = track->GetClusterIndex(i);
3520 Int_t l=(cindex & 0xf0000000) >> 28;
3521 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3523 for (Int_t ind=0;ind<3;ind++){
3525 if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
3527 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3530 Int_t nclusters = track->GetNumberOfClusters();
3531 if (nclusters > 0) //PH Some tracks don't have any cluster
3532 track->SetFakeRatio(double(nwrong)/double(nclusters));
3534 if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
3536 track->SetLabel(tpcLabel);
3540 //------------------------------------------------------------------------
3541 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
3546 //AliITSRecPoint * clist[6];
3547 // Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
3550 track->SetChi2MIP(9,0);
3551 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3552 Int_t cindex = track->GetClusterIndex(i);
3553 Int_t l=(cindex & 0xf0000000) >> 28;
3554 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3555 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3557 for (Int_t ind=0;ind<3;ind++){
3558 if (cl->GetLabel(ind)==lab) isWrong=0;
3560 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3562 //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue; //shared track
3563 //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
3564 //if (l<4&& !(cl->GetType()==1)) continue;
3565 dedx[accepted]= track->GetSampledEdx(i);
3566 //dedx[accepted]= track->fNormQ[l];
3574 TMath::Sort(accepted,dedx,indexes,kFALSE);
3577 Double_t nl=low*accepted, nu =up*accepted;
3579 Float_t sumweight =0;
3580 for (Int_t i=0; i<accepted; i++) {
3582 if (i<nl+0.1) weight = TMath::Max(1.-(nl-i),0.);
3583 if (i>nu-1) weight = TMath::Max(nu-i,0.);
3584 sumamp+= dedx[indexes[i]]*weight;
3587 track->SetdEdx(sumamp/sumweight);
3589 //------------------------------------------------------------------------
3590 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3593 if (fCoefficients) delete []fCoefficients;
3594 fCoefficients = new Float_t[ntracks*48];
3595 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3597 //------------------------------------------------------------------------
3598 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3604 Float_t theta = track->GetTgl();
3605 Float_t phi = track->GetSnp();
3606 phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
3607 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
3608 //printf(" chi2: tr-cl %f %f tr X %f cl X %f\n",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX());
3610 // Take into account the mis-alignment (bring track to cluster plane)
3611 Double_t xTrOrig=track->GetX();
3612 if (!track->PropagateTo(xTrOrig+cluster->GetX(),0.,0.)) return 1000.;
3613 //printf(" chi2: tr-cl %f %f tr X %f cl X %f\n",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX());
3614 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
3615 // Bring the track back to detector plane in ideal geometry
3616 // [mis-alignment will be accounted for in UpdateMI()]
3617 if (!track->PropagateTo(xTrOrig,0.,0.)) return 1000.;
3619 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3620 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3622 chi2+=0.5*TMath::Min(delta/2,2.);
3623 chi2+=2.*cluster->GetDeltaProbability();
3626 track->SetNy(layer,ny);
3627 track->SetNz(layer,nz);
3628 track->SetSigmaY(layer,erry);
3629 track->SetSigmaZ(layer, errz);
3630 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3631 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/(1.- track->GetSnp()*track->GetSnp())));
3635 //------------------------------------------------------------------------
3636 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3641 Int_t layer = (index & 0xf0000000) >> 28;
3642 track->SetClIndex(layer, index);
3643 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3644 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3645 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3646 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3650 if (cl->GetQ()<=0) return 0; // ingore the "virtual" clusters
3653 // Take into account the mis-alignment (bring track to cluster plane)
3654 Double_t xTrOrig=track->GetX();
3655 //Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);printf("gtr %f %f %f\n",trxyz[0],trxyz[1],trxyz[2]);printf("gcl %f %f %f\n",clxyz[0],clxyz[1],clxyz[2]);
3656 //printf(" xtr %f xcl %f\n",track->GetX(),cl->GetX());
3658 if (!track->PropagateTo(xTrOrig+cl->GetX(),0.,0.)) return 0;
3662 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3663 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3666 Int_t updated = track->UpdateMI(&c,chi2,index);
3668 // Bring the track back to detector plane in ideal geometry
3669 if (!track->PropagateTo(xTrOrig,0.,0.)) return 0;
3671 //if(!updated) printf("update failed\n");
3675 //------------------------------------------------------------------------
3676 void AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3679 //DCA sigmas parameterization
3680 //to be paramterized using external parameters in future
3683 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3684 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3686 //------------------------------------------------------------------------
3687 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
3691 Int_t entries = ClusterArray->GetEntriesFast();
3692 if (entries<4) return;
3693 AliITSRecPoint* cluster = (AliITSRecPoint*)ClusterArray->At(0);
3694 Int_t layer = cluster->GetLayer();
3695 if (layer>1) return;
3697 Int_t ncandidates=0;
3698 Float_t r = (layer>0)? 7:4;
3700 for (Int_t i=0;i<entries;i++){
3701 AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(i);
3702 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3703 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3704 index[ncandidates] = i; //candidate to belong to delta electron track
3706 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3707 cl0->SetDeltaProbability(1);
3713 for (Int_t i=0;i<ncandidates;i++){
3714 AliITSRecPoint* cl0 = (AliITSRecPoint*)ClusterArray->At(index[i]);
3715 if (cl0->GetDeltaProbability()>0.8) continue;
3718 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3719 sumy=sumz=sumy2=sumyz=sumw=0.0;
3720 for (Int_t j=0;j<ncandidates;j++){
3722 AliITSRecPoint* cl1 = (AliITSRecPoint*)ClusterArray->At(index[j]);
3724 Float_t dz = cl0->GetZ()-cl1->GetZ();
3725 Float_t dy = cl0->GetY()-cl1->GetY();
3726 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3727 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3728 y[ncl] = cl1->GetY();
3729 z[ncl] = cl1->GetZ();
3730 sumy+= y[ncl]*weight;
3731 sumz+= z[ncl]*weight;
3732 sumy2+=y[ncl]*y[ncl]*weight;
3733 sumyz+=y[ncl]*z[ncl]*weight;
3738 if (ncl<4) continue;
3739 Float_t det = sumw*sumy2 - sumy*sumy;
3741 if (TMath::Abs(det)>0.01){
3742 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3743 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3744 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3747 Float_t z0 = sumyz/sumy;
3748 delta = TMath::Abs(cl0->GetZ()-z0);
3751 cl0->SetDeltaProbability(1-20.*delta);
3755 //------------------------------------------------------------------------
3756 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3760 track->UpdateESDtrack(flags);
3761 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3762 if (oldtrack) delete oldtrack;
3763 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3764 if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3765 printf("Problem\n");
3768 //------------------------------------------------------------------------
3769 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3771 // Get nearest upper layer close to the point xr.
3772 // rough approximation
3774 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3775 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3777 for (Int_t i=0;i<6;i++){
3778 if (radius<kRadiuses[i]){
3785 //------------------------------------------------------------------------
3786 void AliITStrackerMI::UpdateTPCV0(AliESDEvent *event){
3788 //try to update, or reject TPC V0s
3790 Int_t nv0s = event->GetNumberOfV0s();
3791 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3793 for (Int_t i=0;i<nv0s;i++){
3794 AliESDv0 * vertex = event->GetV0(i);
3795 Int_t ip = vertex->GetIndex(0);
3796 Int_t im = vertex->GetIndex(1);
3798 TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
3799 TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
3800 AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
3801 AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
3805 if (trackp->GetNumberOfClusters()+trackp->GetNDeadZone()>5.5){
3806 if (trackp->GetConstrain()&&trackp->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3807 if (!trackp->GetConstrain()&&trackp->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3812 if (trackm->GetNumberOfClusters()+trackm->GetNDeadZone()>5.5){
3813 if (trackm->GetConstrain()&&trackm->GetChi2MIP(0)<3) vertex->SetStatus(-100);
3814 if (!trackm->GetConstrain()&&trackm->GetChi2MIP(0)<2) vertex->SetStatus(-100);
3817 if (vertex->GetStatus()==-100) continue;
3819 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
3820 Int_t clayer = GetNearestLayer(xrp); //I.B.
3821 vertex->SetNBefore(clayer); //
3822 vertex->SetChi2Before(9*clayer); //
3823 vertex->SetNAfter(6-clayer); //
3824 vertex->SetChi2After(0); //
3826 if (clayer >1 ){ // calculate chi2 before vertex
3827 Float_t chi2p = 0, chi2m=0;
3830 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3831 if (trackp->GetClIndex(ilayer)>0){
3832 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3833 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3844 for (Int_t ilayer=0;ilayer<clayer;ilayer++){
3845 if (trackm->GetClIndex(ilayer)>0){
3846 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3847 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3856 vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
3857 if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10); // track exist before vertex
3860 if (clayer < 5 ){ // calculate chi2 after vertex
3861 Float_t chi2p = 0, chi2m=0;
3863 if (trackp&&TMath::Abs(trackp->GetTgl())<1.){
3864 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3865 if (trackp->GetClIndex(ilayer)>0){
3866 chi2p+=trackp->GetDy(ilayer)*trackp->GetDy(ilayer)/(trackp->GetSigmaY(ilayer)*trackp->GetSigmaY(ilayer))+
3867 trackp->GetDz(ilayer)*trackp->GetDz(ilayer)/(trackp->GetSigmaZ(ilayer)*trackp->GetSigmaZ(ilayer));
3877 if (trackm&&TMath::Abs(trackm->GetTgl())<1.){
3878 for (Int_t ilayer=clayer;ilayer<6;ilayer++){
3879 if (trackm->GetClIndex(ilayer)>0){
3880 chi2m+=trackm->GetDy(ilayer)*trackm->GetDy(ilayer)/(trackm->GetSigmaY(ilayer)*trackm->GetSigmaY(ilayer))+
3881 trackm->GetDz(ilayer)*trackm->GetDz(ilayer)/(trackm->GetSigmaZ(ilayer)*trackm->GetSigmaZ(ilayer));
3890 vertex->SetChi2After(TMath::Max(chi2p,chi2m));
3891 if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20); // track not found in ITS
3896 //------------------------------------------------------------------------
3897 void AliITStrackerMI::FindV02(AliESDEvent *event)
3902 // Cuts on DCA - R dependent
3903 // max distance DCA between 2 tracks cut
3904 // maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
3906 const Float_t kMaxDist0 = 0.1;
3907 const Float_t kMaxDist1 = 0.1;
3908 const Float_t kMaxDist = 1;
3909 const Float_t kMinPointAngle = 0.85;
3910 const Float_t kMinPointAngle2 = 0.99;
3911 const Float_t kMinR = 0.5;
3912 const Float_t kMaxR = 220;
3913 //const Float_t kCausality0Cut = 0.19;
3914 //const Float_t kLikelihood01Cut = 0.25;
3915 //const Float_t kPointAngleCut = 0.9996;
3916 const Float_t kCausality0Cut = 0.19;
3917 const Float_t kLikelihood01Cut = 0.45;
3918 const Float_t kLikelihood1Cut = 0.5;
3919 const Float_t kCombinedCut = 0.55;
3923 TTreeSRedirector &cstream = *fDebugStreamer;
3924 Int_t ntracks = event->GetNumberOfTracks();
3925 Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
3926 fOriginal.Expand(ntracks);
3927 fTrackHypothesys.Expand(ntracks);
3928 fBestHypothesys.Expand(ntracks);
3930 AliHelix * helixes = new AliHelix[ntracks+2];
3931 TObjArray trackarray(ntracks+2); //array with tracks - with vertex constrain
3932 TObjArray trackarrayc(ntracks+2); //array of "best tracks" - without vertex constrain
3933 TObjArray trackarrayl(ntracks+2); //array of "longest tracks" - without vertex constrain
3934 Bool_t * forbidden = new Bool_t [ntracks+2];
3935 Int_t *itsmap = new Int_t [ntracks+2];
3936 Float_t *dist = new Float_t[ntracks+2];
3937 Float_t *normdist0 = new Float_t[ntracks+2];
3938 Float_t *normdist1 = new Float_t[ntracks+2];
3939 Float_t *normdist = new Float_t[ntracks+2];
3940 Float_t *norm = new Float_t[ntracks+2];
3941 Float_t *maxr = new Float_t[ntracks+2];
3942 Float_t *minr = new Float_t[ntracks+2];
3943 Float_t *minPointAngle= new Float_t[ntracks+2];
3945 AliV0 *pvertex = new AliV0;
3946 AliITStrackMI * dummy= new AliITStrackMI;
3948 AliITStrackMI trackat0; //temporary track for DCA calculation
3950 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3952 // make ITS - ESD map
3954 for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
3955 itsmap[itrack] = -1;
3956 forbidden[itrack] = kFALSE;
3957 maxr[itrack] = kMaxR;
3958 minr[itrack] = kMinR;
3959 minPointAngle[itrack] = kMinPointAngle;
3961 for (Int_t itrack=0;itrack<nitstracks;itrack++){
3962 AliITStrackMI * original = (AliITStrackMI*)(fOriginal.At(itrack));
3963 Int_t esdindex = original->GetESDtrack()->GetID();
3964 itsmap[esdindex] = itrack;
3967 // create ITS tracks from ESD tracks if not done before
3969 for (Int_t itrack=0;itrack<ntracks;itrack++){
3970 if (itsmap[itrack]>=0) continue;
3971 AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
3972 //tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
3973 //tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ();
3974 tpctrack->GetDZ(GetX(),GetY(),GetZ(),tpctrack->GetDP()); //I.B.
3975 if (tpctrack->GetD(0)<20 && tpctrack->GetD(1)<20){
3976 // tracks which can reach inner part of ITS
3977 // propagate track to outer its volume - with correction for material
3978 CorrectForTPCtoITSDeadZoneMaterial(tpctrack);
3980 itsmap[itrack] = nitstracks;
3981 fOriginal.AddAt(tpctrack,nitstracks);
3985 // fill temporary arrays
3987 for (Int_t itrack=0;itrack<ntracks;itrack++){
3988 AliESDtrack * esdtrack = event->GetTrack(itrack);
3989 Int_t itsindex = itsmap[itrack];
3990 AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
3991 if (!original) continue;
3992 AliITStrackMI *bestConst = 0;
3993 AliITStrackMI *bestLong = 0;
3994 AliITStrackMI *best = 0;
3997 TObjArray * array = (TObjArray*) fTrackHypothesys.At(itsindex);
3998 Int_t hentries = (array==0) ? 0 : array->GetEntriesFast();
3999 // Get best track with vertex constrain
4000 for (Int_t ih=0;ih<hentries;ih++){
4001 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4002 if (!trackh->GetConstrain()) continue;
4003 if (!bestConst) bestConst = trackh;
4004 if (trackh->GetNumberOfClusters()>5.0){
4005 bestConst = trackh; // full track - with minimal chi2
4008 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()) continue;
4012 // Get best long track without vertex constrain and best track without vertex constrain
4013 for (Int_t ih=0;ih<hentries;ih++){
4014 AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
4015 if (trackh->GetConstrain()) continue;
4016 if (!best) best = trackh;
4017 if (!bestLong) bestLong = trackh;
4018 if (trackh->GetNumberOfClusters()>5.0){
4019 bestLong = trackh; // full track - with minimal chi2
4022 if (trackh->GetNumberOfClusters()+trackh->GetNDeadZone()<=bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()) continue;
4027 bestLong = original;
4029 //I.B. trackat0 = *bestLong;
4030 new (&trackat0) AliITStrackMI(*bestLong);
4031 Double_t xx,yy,zz,alpha;
4032 bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
4033 alpha = TMath::ATan2(yy,xx);
4034 trackat0.Propagate(alpha,0);
4035 // calculate normalized distances to the vertex
4037 Float_t ptfac = (1.+100.*TMath::Abs(trackat0.GetC()));
4038 if ( bestLong->GetNumberOfClusters()>3 ){
4039 dist[itsindex] = trackat0.GetY();
4040 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4041 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4042 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4043 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4045 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<6) normdist[itsindex]*=2.;
4046 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<5) normdist[itsindex]*=2.;
4047 if (bestLong->GetNumberOfClusters()+bestLong->GetNDeadZone()<4) normdist[itsindex]*=2.;
4049 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<6) normdist[itsindex]*=1.5;
4050 if (bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()<5) normdist[itsindex]*=1.5;
4054 if (bestConst&&bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5){
4055 dist[itsindex] = bestConst->GetD(0);
4056 norm[itsindex] = bestConst->GetDnorm(0);
4057 normdist0[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4058 normdist1[itsindex] = TMath::Abs(bestConst->GetD(0)/norm[itsindex]);
4059 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4061 dist[itsindex] = trackat0.GetY();
4062 norm[itsindex] = ptfac*TMath::Sqrt(trackat0.GetSigmaY2());
4063 normdist0[itsindex] = TMath::Abs(trackat0.GetY()/norm[itsindex]);
4064 normdist1[itsindex] = TMath::Abs((trackat0.GetZ()-primvertex[2])/(ptfac*TMath::Sqrt(trackat0.GetSigmaZ2())));
4065 normdist[itsindex] = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
4066 if (TMath::Abs(trackat0.GetTgl())>1.05){
4067 if (normdist[itsindex]<3) forbidden[itsindex]=kTRUE;
4068 if (normdist[itsindex]>3) {
4069 minr[itsindex] = TMath::Max(Float_t(40.),minr[itsindex]);
4075 //-----------------------------------------------------------
4076 //Forbid primary track candidates -
4078 //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
4079 //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
4080 //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
4081 //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
4082 //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
4083 //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
4084 //-----------------------------------------------------------
4086 if (bestLong->GetNumberOfClusters()<4 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>4.5) forbidden[itsindex]=kTRUE;
4087 if (normdist[itsindex]<3 && bestConst->GetNumberOfClusters()+bestConst->GetNDeadZone()>5.5) forbidden[itsindex]=kTRUE;
4088 if (normdist[itsindex]<2 && bestConst->GetClIndex(0)>0 && bestConst->GetClIndex(1)>0 ) forbidden[itsindex]=kTRUE;
4089 if (normdist[itsindex]<1 && bestConst->GetClIndex(0)>0) forbidden[itsindex]=kTRUE;
4090 if (normdist[itsindex]<4 && bestConst->GetNormChi2(0)<2) forbidden[itsindex]=kTRUE;
4091 if (normdist[itsindex]<5 && bestConst->GetNormChi2(0)<1) forbidden[itsindex]=kTRUE;
4092 if (bestConst->GetNormChi2(0)<2.5) {
4093 minPointAngle[itsindex]= 0.9999;
4094 maxr[itsindex] = 10;
4098 //forbid daughter kink candidates
4100 if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
4101 Bool_t isElectron = kTRUE;
4102 Bool_t isProton = kTRUE;
4104 esdtrack->GetESDpid(pid);
4105 for (Int_t i=1;i<5;i++){
4106 if (pid[0]<pid[i]) isElectron= kFALSE;
4107 if (pid[4]<pid[i]) isProton= kFALSE;
4110 forbidden[itsindex]=kFALSE;
4111 normdist[itsindex]*=-1;
4114 if (normdist[itsindex]>2) forbidden[itsindex]=kFALSE;
4115 normdist[itsindex]*=-1;
4119 // Causality cuts in TPC volume
4121 if (esdtrack->GetTPCdensity(0,10) >0.6) maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
4122 if (esdtrack->GetTPCdensity(10,30)>0.6) maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
4123 if (esdtrack->GetTPCdensity(20,40)>0.6) maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
4124 if (esdtrack->GetTPCdensity(30,50)>0.6) maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
4126 if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->GetNumberOfClusters()<3) minr[itsindex]=100;
4132 "Tr1.="<<((bestConst)? bestConst:dummy)<<
4134 "Tr3.="<<&trackat0<<
4136 "Dist="<<dist[itsindex]<<
4137 "ND0="<<normdist0[itsindex]<<
4138 "ND1="<<normdist1[itsindex]<<
4139 "ND="<<normdist[itsindex]<<
4140 "Pz="<<primvertex[2]<<
4141 "Forbid="<<forbidden[itsindex]<<
4145 trackarray.AddAt(best,itsindex);
4146 trackarrayc.AddAt(bestConst,itsindex);
4147 trackarrayl.AddAt(bestLong,itsindex);
4148 new (&helixes[itsindex]) AliHelix(*best);
4153 // first iterration of V0 finder
4155 for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
4156 Int_t itrack0 = itsmap[iesd0];
4157 if (forbidden[itrack0]) continue;
4158 AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
4159 if (!btrack0) continue;
4160 if (btrack0->GetSign()>0) continue;
4161 AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
4163 for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
4164 Int_t itrack1 = itsmap[iesd1];
4165 if (forbidden[itrack1]) continue;
4167 AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1);
4168 if (!btrack1) continue;
4169 if (btrack1->GetSign()<0) continue;
4170 Bool_t isGold = kFALSE;
4171 if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
4174 AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
4175 AliHelix &h1 = helixes[itrack0];
4176 AliHelix &h2 = helixes[itrack1];
4178 // find linear distance
4183 Double_t phase[2][2],radius[2];
4184 Int_t points = h1.GetRPHIintersections(h2, phase, radius);
4185 if (points==0) continue;
4186 Double_t delta[2]={1000000,1000000};
4188 h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
4190 if (radius[1]<rmin) rmin = radius[1];
4191 h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
4193 rmin = TMath::Sqrt(rmin);
4194 Double_t distance = 0;
4195 Double_t radiusC = 0;
4197 if (points==1 || delta[0]<delta[1]){
4198 distance = TMath::Sqrt(delta[0]);
4199 radiusC = TMath::Sqrt(radius[0]);
4201 distance = TMath::Sqrt(delta[1]);
4202 radiusC = TMath::Sqrt(radius[1]);
4205 if (radiusC<TMath::Max(minr[itrack0],minr[itrack1])) continue;
4206 if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1])) continue;
4207 Float_t maxDist = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));
4208 if (distance>maxDist) continue;
4209 Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
4210 if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
4213 // Double_t distance = TestV0(h1,h2,pvertex,rmin);
4215 // if (distance>maxDist) continue;
4216 // if (pvertex->GetRr()<kMinR) continue;
4217 // if (pvertex->GetRr()>kMaxR) continue;
4218 AliITStrackMI * track0=btrack0;
4219 AliITStrackMI * track1=btrack1;
4220 // if (pvertex->GetRr()<3.5){
4222 //use longest tracks inside the pipe
4223 track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
4224 track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
4228 pvertex->SetParamN(*track0);
4229 pvertex->SetParamP(*track1);
4230 pvertex->Update(primvertex);
4231 pvertex->SetClusters(track0->ClIndex(),track1->ClIndex()); // register clusters
4233 if (pvertex->GetRr()<kMinR) continue;
4234 if (pvertex->GetRr()>kMaxR) continue;
4235 if (pvertex->GetV0CosineOfPointingAngle()<kMinPointAngle) continue;
4236 //Bo: if (pvertex->GetDist2()>maxDist) continue;
4237 if (pvertex->GetDcaV0Daughters()>maxDist) continue;
4238 //Bo: pvertex->SetLab(0,track0->GetLabel());
4239 //Bo: pvertex->SetLab(1,track1->GetLabel());
4240 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4241 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4243 AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
4244 AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
4248 TObjArray * array0b = (TObjArray*)fBestHypothesys.At(itrack0);
4249 if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->GetTgl())<1.1) {
4250 fCurrentEsdTrack = itrack0;
4251 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
4253 TObjArray * array1b = (TObjArray*)fBestHypothesys.At(itrack1);
4254 if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->GetTgl())<1.1) {
4255 fCurrentEsdTrack = itrack1;
4256 FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
4259 AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);
4260 AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
4261 AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);
4262 AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
4264 Float_t minchi2before0=16;
4265 Float_t minchi2before1=16;
4266 Float_t minchi2after0 =16;
4267 Float_t minchi2after1 =16;
4268 Double_t xrp[3]; pvertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
4269 Int_t maxLayer = GetNearestLayer(xrp); //I.B.
4271 if (array0b) for (Int_t i=0;i<5;i++){
4272 // best track after vertex
4273 AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
4274 if (!btrack) continue;
4275 if (btrack->GetNumberOfClusters()>track0l->GetNumberOfClusters()) track0l = btrack;
4276 // if (btrack->fX<pvertex->GetRr()-2.-0.5/(0.1+pvertex->GetAnglep()[2])) {
4277 if (btrack->GetX()<pvertex->GetRr()-2.) {
4278 if ( (maxLayer>i+2|| (i==0)) && btrack->GetNumberOfClusters()==(6-i)&&i<3){
4281 if (maxLayer<3){ // take prim vertex as additional measurement
4282 if (normdist[itrack0]>0 && htrackc0){
4283 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack0]*normdist[itrack0],16.);
4285 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack0]*normdist[itrack0]+3.),16.);
4289 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4291 if (!btrack->GetClIndex(ilayer)){
4295 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4296 for (Int_t itrack=0;itrack<4;itrack++){
4297 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack0){
4298 sumchi2+=18.; //shared cluster
4302 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4303 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4307 if (sumchi2<minchi2before0) minchi2before0=sumchi2;
4309 continue; //safety space - Geo manager will give exact layer
4312 minchi2after0 = btrack->GetNormChi2(i);
4315 if (array1b) for (Int_t i=0;i<5;i++){
4316 // best track after vertex
4317 AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
4318 if (!btrack) continue;
4319 if (btrack->GetNumberOfClusters()>track1l->GetNumberOfClusters()) track1l = btrack;
4320 // if (btrack->fX<pvertex->GetRr()-2-0.5/(0.1+pvertex->GetAnglep()[2])){
4321 if (btrack->GetX()<pvertex->GetRr()-2){
4322 if ((maxLayer>i+2 || (i==0))&&btrack->GetNumberOfClusters()==(6-i)&&(i<3)){
4325 if (maxLayer<3){ // take prim vertex as additional measurement
4326 if (normdist[itrack1]>0 && htrackc1){
4327 sumchi2 += TMath::Min((3.-maxLayer)*normdist[itrack1]*normdist[itrack1],16.);
4329 sumchi2 += TMath::Min((3.-maxLayer)*(3*normdist[itrack1]*normdist[itrack1]+3.),16.);
4333 for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
4335 if (!btrack->GetClIndex(ilayer)){
4339 Int_t c=( btrack->GetClIndex(ilayer) & 0x0fffffff);
4340 for (Int_t itrack=0;itrack<4;itrack++){
4341 if (fgLayers[ilayer].GetClusterTracks(itrack,c)>=0 && fgLayers[ilayer].GetClusterTracks(itrack,c)!=itrack1){
4342 sumchi2+=18.; //shared cluster
4346 sumchi2+=btrack->GetDy(ilayer)*btrack->GetDy(ilayer)/(btrack->GetSigmaY(ilayer)*btrack->GetSigmaY(ilayer));
4347 sumchi2+=btrack->GetDz(ilayer)*btrack->GetDz(ilayer)/(btrack->GetSigmaZ(ilayer)*btrack->GetSigmaZ(ilayer));
4351 if (sumchi2<minchi2before1) minchi2before1=sumchi2;
4353 continue; //safety space - Geo manager will give exact layer
4356 minchi2after1 = btrack->GetNormChi2(i);
4360 // position resolution - used for DCA cut
4361 Float_t sigmad = track0b->GetSigmaY2()+track0b->GetSigmaZ2()+track1b->GetSigmaY2()+track1b->GetSigmaZ2()+
4362 (track0b->GetX()-pvertex->GetRr())*(track0b->GetX()-pvertex->GetRr())*(track0b->GetSigmaSnp2()+track0b->GetSigmaTgl2())+
4363 (track1b->GetX()-pvertex->GetRr())*(track1b->GetX()-pvertex->GetRr())*(track1b->GetSigmaSnp2()+track1b->GetSigmaTgl2());
4364 sigmad =TMath::Sqrt(sigmad)+0.04;
4365 if (pvertex->GetRr()>50){
4366 Double_t cov0[15],cov1[15];
4367 track0b->GetESDtrack()->GetInnerExternalCovariance(cov0);
4368 track1b->GetESDtrack()->GetInnerExternalCovariance(cov1);
4369 sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
4370 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
4371 (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
4372 sigmad =TMath::Sqrt(sigmad)+0.05;
4376 vertex2.SetParamN(*track0b);
4377 vertex2.SetParamP(*track1b);
4378 vertex2.Update(primvertex);
4379 //Bo: if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4380 if (vertex2.GetDcaV0Daughters()<=pvertex->GetDcaV0Daughters()&&(vertex2.GetV0CosineOfPointingAngle()>=pvertex->GetV0CosineOfPointingAngle())){
4381 pvertex->SetParamN(*track0b);
4382 pvertex->SetParamP(*track1b);
4383 pvertex->Update(primvertex);
4384 pvertex->SetClusters(track0b->ClIndex(),track1b->ClIndex()); // register clusters
4385 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
4386 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
4388 pvertex->SetDistSigma(sigmad);
4389 //Bo: pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);
4390 pvertex->SetNormDCAPrim(normdist[itrack0],normdist[itrack1]);
4392 // define likelihhod and causalities
4394 Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;
4396 Float_t fnorm0 = normdist[itrack0];
4397 if (fnorm0<0) fnorm0*=-3;
4398 Float_t fnorm1 = normdist[itrack1];
4399 if (fnorm1<0) fnorm1*=-3;
4400 if (pvertex->GetAnglep()[2]>0.1 || (pvertex->GetRr()<10.5)&& pvertex->GetAnglep()[2]>0.05 || pvertex->GetRr()<3){
4401 pb0 = TMath::Exp(-TMath::Min(fnorm0,Float_t(16.))/12.);
4402 pb1 = TMath::Exp(-TMath::Min(fnorm1,Float_t(16.))/12.);
4404 pvertex->SetChi2Before(normdist[itrack0]);
4405 pvertex->SetChi2After(normdist[itrack1]);
4406 pvertex->SetNAfter(0);
4407 pvertex->SetNBefore(0);
4409 pvertex->SetChi2Before(minchi2before0);
4410 pvertex->SetChi2After(minchi2before1);
4411 if (pvertex->GetAnglep()[2]>0.1 || ( pvertex->GetRr()<10.5 && pvertex->GetAnglep()[2]>0.05) || pvertex->GetRr()<3){
4412 pb0 = TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
4413 pb1 = TMath::Exp(-TMath::Min(minchi2before1,Float_t(16))/12.);
4415 pvertex->SetNAfter(maxLayer);
4416 pvertex->SetNBefore(maxLayer);
4418 if (pvertex->GetRr()<90){
4419 pa0 *= TMath::Min(track0->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4420 pa1 *= TMath::Min(track1->GetESDtrack()->GetTPCdensity(0,60),Double_t(1.));
4422 if (pvertex->GetRr()<20){
4423 pa0 *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
4424 pa1 *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
4427 pvertex->SetCausality(pb0,pb1,pa0,pa1);
4429 // Likelihood calculations - apply cuts
4431 Bool_t v0OK = kTRUE;
4432 Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
4433 p12 += pvertex->GetParamN()->GetParameter()[4]*pvertex->GetParamN()->GetParameter()[4];
4434 p12 = TMath::Sqrt(p12); // "mean" momenta
4435 Float_t sigmap0 = 0.0001+0.001/(0.1+pvertex->GetRr());
4436 Float_t sigmap = 0.5*sigmap0*(0.6+0.4*p12); // "resolution: of point angle - as a function of radius and momenta
4438 Float_t causalityA = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
4439 Float_t causalityB = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Double_t(0.7))*
4440 TMath::Min(pvertex->GetCausalityP()[3],Double_t(0.7)));
4442 //Bo: Float_t likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
4443 Float_t lDistNorm = pvertex->GetDcaV0Daughters()/pvertex->GetDistSigma();
4444 Float_t likelihood0 = (TMath::Exp(-lDistNorm)+0.1) *(pvertex->GetDcaV0Daughters()<0.5)*(lDistNorm<5);
4446 Float_t likelihood1 = TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/sigmap)+
4447 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(4.*sigmap))+
4448 0.4*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/(8.*sigmap))+
4449 0.1*TMath::Exp(-(1.0001-pvertex->GetV0CosineOfPointingAngle())/0.01);
4451 if (causalityA<kCausality0Cut) v0OK = kFALSE;
4452 if (TMath::Sqrt(likelihood0*likelihood1)<kLikelihood01Cut) v0OK = kFALSE;
4453 if (likelihood1<kLikelihood1Cut) v0OK = kFALSE;
4454 if (TMath::Power(likelihood0*likelihood1*causalityB,0.33)<kCombinedCut) v0OK = kFALSE;
4459 Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
4461 "Tr0.="<<track0<< //best without constrain
4462 "Tr1.="<<track1<< //best without constrain
4463 "Tr0B.="<<track0b<< //best without constrain after vertex
4464 "Tr1B.="<<track1b<< //best without constrain after vertex
4465 "Tr0C.="<<htrackc0<< //best with constrain if exist
4466 "Tr1C.="<<htrackc1<< //best with constrain if exist
4467 "Tr0L.="<<track0l<< //longest best
4468 "Tr1L.="<<track1l<< //longest best
4469 "Esd0.="<<track0->GetESDtrack()<< // esd track0 params
4470 "Esd1.="<<track1->GetESDtrack()<< // esd track1 params
4471 "V0.="<<pvertex<< //vertex properties
4472 "V0b.="<<&vertex2<< //vertex properties at "best" track
4473 "ND0="<<normdist[itrack0]<< //normalize distance for track0
4474 "ND1="<<normdist[itrack1]<< //normalize distance for track1
4476 // "RejectBase="<<rejectBase<< //rejection in First itteration
4482 //if (rejectBase) continue;
4484 pvertex->SetStatus(0);
4485 // if (rejectBase) {
4486 // pvertex->SetStatus(-100);
4488 if (pvertex->GetV0CosineOfPointingAngle()>kMinPointAngle2) {
4489 //Bo: pvertex->SetESDindexes(track0->GetESDtrack()->GetID(),track1->GetESDtrack()->GetID());
4490 pvertex->SetIndex(0,track0->GetESDtrack()->GetID());//Bo: consistency 0 for neg
4491 pvertex->SetIndex(1,track1->GetESDtrack()->GetID());//Bo: consistency 1 for pos
4493 // AliV0vertex vertexjuri(*track0,*track1);
4494 // vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
4495 // event->AddV0(&vertexjuri);
4496 pvertex->SetStatus(100);
4498 pvertex->SetOnFlyStatus(kTRUE);
4499 pvertex->ChangeMassHypothesis(kK0Short);
4500 event->AddV0(pvertex);
4506 // delete temporary arrays
4509 delete[] minPointAngle;
4521 //------------------------------------------------------------------------
4522 void AliITStrackerMI::RefitV02(AliESDEvent *event)
4525 //try to refit V0s in the third path of the reconstruction
4527 TTreeSRedirector &cstream = *fDebugStreamer;
4529 Int_t nv0s = event->GetNumberOfV0s();
4530 Float_t primvertex[3]={GetX(),GetY(),GetZ()};
4532 for (Int_t iv0 = 0; iv0<nv0s;iv0++){
4533 AliV0 * v0mi = (AliV0*)event->GetV0(iv0);
4534 if (!v0mi) continue;
4535 Int_t itrack0 = v0mi->GetIndex(0);
4536 Int_t itrack1 = v0mi->GetIndex(1);
4537 AliESDtrack *esd0 = event->GetTrack(itrack0);
4538 AliESDtrack *esd1 = event->GetTrack(itrack1);
4539 if (!esd0||!esd1) continue;
4540 AliITStrackMI tpc0(*esd0);
4541 AliITStrackMI tpc1(*esd1);
4542 Double_t x,y,z; v0mi->GetXYZ(x,y,z); //I.B.
4543 Double_t alpha =TMath::ATan2(y,x); //I.B.
4544 if (v0mi->GetRr()>85){
4545 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4546 v0temp.SetParamN(tpc0);
4547 v0temp.SetParamP(tpc1);
4548 v0temp.Update(primvertex);
4549 if (kFALSE) cstream<<"Refit"<<
4551 "V0refit.="<<&v0temp<<
4555 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4556 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4557 v0mi->SetParamN(tpc0);
4558 v0mi->SetParamP(tpc1);
4559 v0mi->Update(primvertex);
4564 if (v0mi->GetRr()>35){
4565 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4566 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4567 if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4568 v0temp.SetParamN(tpc0);
4569 v0temp.SetParamP(tpc1);
4570 v0temp.Update(primvertex);
4571 if (kFALSE) cstream<<"Refit"<<
4573 "V0refit.="<<&v0temp<<
4577 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4578 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4579 v0mi->SetParamN(tpc0);
4580 v0mi->SetParamP(tpc1);
4581 v0mi->Update(primvertex);
4586 CorrectForTPCtoITSDeadZoneMaterial(&tpc0);
4587 CorrectForTPCtoITSDeadZoneMaterial(&tpc1);
4588 // if (tpc0.Propagate(alpha,v0mi->GetRr())&&tpc1.Propagate(alpha,v0mi->GetRr())){
4589 if (RefitAt(v0mi->GetRr(),&tpc0, v0mi->GetClusters(0)) && RefitAt(v0mi->GetRr(),&tpc1, v0mi->GetClusters(1))){
4590 v0temp.SetParamN(tpc0);
4591 v0temp.SetParamP(tpc1);
4592 v0temp.Update(primvertex);
4593 if (kFALSE) cstream<<"Refit"<<
4595 "V0refit.="<<&v0temp<<
4599 //Bo: if (v0temp.GetDist2()<v0mi->GetDist2() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4600 if (v0temp.GetDcaV0Daughters()<v0mi->GetDcaV0Daughters() || v0temp.GetV0CosineOfPointingAngle()>v0mi->GetV0CosineOfPointingAngle()){
4601 v0mi->SetParamN(tpc0);
4602 v0mi->SetParamP(tpc1);
4603 v0mi->Update(primvertex);
4608 //------------------------------------------------------------------------
4609 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4610 //--------------------------------------------------------------------
4611 // Fill a look-up table with mean material
4612 //--------------------------------------------------------------------
4616 Double_t point1[3],point2[3];
4617 Double_t phi,cosphi,sinphi,z;
4618 // 0-5 layers, 6 pipe, 7-8 shields
4619 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4620 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4622 Int_t ifirst=0,ilast=0;
4623 if(material.Contains("Pipe")) {
4625 } else if(material.Contains("Shields")) {
4627 } else if(material.Contains("Layers")) {
4630 Error("BuildMaterialLUT","Wrong layer name\n");
4633 for(Int_t imat=ifirst; imat<=ilast; imat++) {
4634 Double_t param[5]={0.,0.,0.,0.,0.};
4635 for (Int_t i=0; i<n; i++) {
4636 phi = 2.*TMath::Pi()*gRandom->Rndm();
4637 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4638 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4639 point1[0] = rmin[imat]*cosphi;
4640 point1[1] = rmin[imat]*sinphi;
4642 point2[0] = rmax[imat]*cosphi;
4643 point2[1] = rmax[imat]*sinphi;
4645 AliTracker::MeanMaterialBudget(point1,point2,mparam);
4646 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4648 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4650 fxOverX0Layer[imat] = param[1];
4651 fxTimesRhoLayer[imat] = param[0]*param[4];
4652 } else if(imat==6) {
4653 fxOverX0Pipe = param[1];
4654 fxTimesRhoPipe = param[0]*param[4];
4655 } else if(imat==7) {
4656 fxOverX0Shield[0] = param[1];
4657 fxTimesRhoShield[0] = param[0]*param[4];
4658 } else if(imat==8) {
4659 fxOverX0Shield[1] = param[1];
4660 fxTimesRhoShield[1] = param[0]*param[4];
4664 printf("%s\n",material.Data());
4665 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4666 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4667 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4668 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4669 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4670 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4671 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4672 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4673 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4677 //------------------------------------------------------------------------
4678 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4679 TString direction) {
4680 //-------------------------------------------------------------------
4681 // Propagate beyond beam pipe and correct for material
4682 // (material budget in different ways according to fUseTGeo value)
4683 //-------------------------------------------------------------------
4685 // Define budget mode:
4686 // 0: material from AliITSRecoParam (hard coded)
4687 // 1: material from TGeo (on the fly)
4688 // 2: material from lut
4689 // 3: material from TGeo (same for all hypotheses)
4702 if(fTrackingPhase.Contains("Clusters2Tracks"))
4703 { mode=3; } else { mode=1; }
4706 if(fTrackingPhase.Contains("Clusters2Tracks"))
4707 { mode=3; } else { mode=2; }
4713 if(fTrackingPhase.Contains("Default")) mode=0;
4715 Int_t index=fCurrentEsdTrack;
4717 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4718 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4720 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4722 Double_t xOverX0,x0,lengthTimesMeanDensity;
4723 Bool_t anglecorr=kTRUE;
4727 xOverX0 = AliITSRecoParam::GetdPipe();
4728 x0 = AliITSRecoParam::GetX0Be();
4729 lengthTimesMeanDensity = xOverX0*x0;
4732 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4736 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4737 xOverX0 = fxOverX0Pipe;
4738 lengthTimesMeanDensity = fxTimesRhoPipe;
4741 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4742 if(fxOverX0PipeTrks[index]<0) {
4743 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4744 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4745 (1.-t->GetSnp()*t->GetSnp()));
4746 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4747 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4750 xOverX0 = fxOverX0PipeTrks[index];
4751 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4755 lengthTimesMeanDensity *= dir;
4757 if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4758 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4762 //------------------------------------------------------------------------
4763 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4765 TString direction) {
4766 //-------------------------------------------------------------------
4767 // Propagate beyond SPD or SDD shield and correct for material
4768 // (material budget in different ways according to fUseTGeo value)
4769 //-------------------------------------------------------------------
4771 // Define budget mode:
4772 // 0: material from AliITSRecoParam (hard coded)
4773 // 1: material from TGeo (on the fly)
4774 // 2: material from lut
4775 // 3: material from TGeo (same for all hypotheses)
4788 if(fTrackingPhase.Contains("Clusters2Tracks"))
4789 { mode=3; } else { mode=1; }
4792 if(fTrackingPhase.Contains("Clusters2Tracks"))
4793 { mode=3; } else { mode=2; }
4799 if(fTrackingPhase.Contains("Default")) mode=0;
4801 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4803 Int_t shieldindex=0;
4804 if (shield.Contains("SDD")) { // SDDouter
4805 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4807 } else if (shield.Contains("SPD")) { // SPDouter
4808 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4811 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4815 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4817 Int_t index=2*fCurrentEsdTrack+shieldindex;
4819 Double_t xOverX0,x0,lengthTimesMeanDensity;
4820 Bool_t anglecorr=kTRUE;
4824 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4825 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4826 lengthTimesMeanDensity = xOverX0*x0;
4829 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4833 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4834 xOverX0 = fxOverX0Shield[shieldindex];
4835 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4838 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4839 if(fxOverX0ShieldTrks[index]<0) {
4840 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4841 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4842 (1.-t->GetSnp()*t->GetSnp()));
4843 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4844 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4847 xOverX0 = fxOverX0ShieldTrks[index];
4848 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4852 lengthTimesMeanDensity *= dir;
4854 if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4855 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4859 //------------------------------------------------------------------------
4860 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4862 Double_t oldGlobXYZ[3],
4863 TString direction) {
4864 //-------------------------------------------------------------------
4865 // Propagate beyond layer and correct for material
4866 // (material budget in different ways according to fUseTGeo value)
4867 //-------------------------------------------------------------------
4869 // Define budget mode:
4870 // 0: material from AliITSRecoParam (hard coded)
4871 // 1: material from TGeo (on the fly)
4872 // 2: material from lut
4873 // 3: material from TGeo (same for all hypotheses)
4886 if(fTrackingPhase.Contains("Clusters2Tracks"))
4887 { mode=3; } else { mode=1; }
4890 if(fTrackingPhase.Contains("Clusters2Tracks"))
4891 { mode=3; } else { mode=2; }
4897 if(fTrackingPhase.Contains("Default")) mode=0;
4899 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4901 Double_t r=fgLayers[layerindex].GetR();
4902 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4904 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4906 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4908 Int_t index=6*fCurrentEsdTrack+layerindex;
4910 // Bring the track beyond the material
4911 if (!t->AliExternalTrackParam::PropagateTo(xToGo,GetBz())) return 0;
4912 Double_t globXYZ[3];
4915 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4917 Bool_t anglecorr=kTRUE;
4921 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4922 lengthTimesMeanDensity = xOverX0*x0;
4925 AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4926 if(mparam[1]>900000) return 0;
4928 lengthTimesMeanDensity=mparam[0]*mparam[4];
4932 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4933 xOverX0 = fxOverX0Layer[layerindex];
4934 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4937 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4938 if(fxOverX0LayerTrks[index]<0) {
4939 AliTracker::MeanMaterialBudget(oldGlobXYZ,globXYZ,mparam);
4940 if(mparam[1]>900000) return 0;
4941 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4942 (1.-t->GetSnp()*t->GetSnp()));
4943 xOverX0=mparam[1]/angle;
4944 lengthTimesMeanDensity=mparam[0]*mparam[4]/angle;
4945 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0);
4946 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity);
4948 xOverX0 = fxOverX0LayerTrks[index];
4949 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4953 lengthTimesMeanDensity *= dir;
4955 if (!t->CorrectForMeanMaterial(xOverX0,lengthTimesMeanDensity,anglecorr)) return 0;
4959 //------------------------------------------------------------------------
4960 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4961 //-----------------------------------------------------------------
4962 // Initialize LUT for storing material for each prolonged track
4963 //-----------------------------------------------------------------
4964 fxOverX0PipeTrks = new Float_t[ntracks];
4965 fxTimesRhoPipeTrks = new Float_t[ntracks];
4966 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4967 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4968 fxOverX0LayerTrks = new Float_t[ntracks*6];
4969 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4971 for(Int_t i=0; i<ntracks; i++) {
4972 fxOverX0PipeTrks[i] = -1.;
4973 fxTimesRhoPipeTrks[i] = -1.;
4975 for(Int_t j=0; j<ntracks*2; j++) {
4976 fxOverX0ShieldTrks[j] = -1.;
4977 fxTimesRhoShieldTrks[j] = -1.;
4979 for(Int_t k=0; k<ntracks*6; k++) {
4980 fxOverX0LayerTrks[k] = -1.;
4981 fxTimesRhoLayerTrks[k] = -1.;
4988 //------------------------------------------------------------------------
4989 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4990 //-----------------------------------------------------------------
4991 // Delete LUT for storing material for each prolonged track
4992 //-----------------------------------------------------------------
4993 if(fxOverX0PipeTrks) {
4994 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4996 if(fxOverX0ShieldTrks) {
4997 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
5000 if(fxOverX0LayerTrks) {
5001 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
5003 if(fxTimesRhoPipeTrks) {
5004 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
5006 if(fxTimesRhoShieldTrks) {
5007 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
5009 if(fxTimesRhoLayerTrks) {
5010 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
5014 //------------------------------------------------------------------------
5015 Int_t AliITStrackerMI::CheckSkipLayer(AliITStrackMI *track,
5016 Int_t ilayer,Int_t idet) const {
5017 //-----------------------------------------------------------------
5018 // This method is used to decide whether to allow a prolongation
5019 // without clusters, because we want to skip the layer.
5020 // In this case the return value is > 0:
5021 // return 1: the user requested to skip a layer
5022 // return 2: track outside z acceptance of SSD/SDD and will cross both SPD
5023 //-----------------------------------------------------------------
5025 if (AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilayer)) return 1;
5027 if (idet<0 && ilayer>1 && AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
5028 // check if track will cross SPD outer layer
5029 Double_t phiAtSPD2,zAtSPD2;
5030 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
5031 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
5037 //------------------------------------------------------------------------
5038 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
5039 Int_t ilayer,Int_t idet,
5040 Double_t dz,Double_t dy,
5041 Bool_t noClusters) const {
5042 //-----------------------------------------------------------------
5043 // This method is used to decide whether to allow a prolongation
5044 // without clusters, because there is a dead zone in the road.
5045 // In this case the return value is > 0:
5046 // return 1: dead zone at z=0,+-7cm in SPD
5047 // return 2: all road is "bad" (dead or noisy) from the OCDB
5048 // return 3: something "bad" (dead or noisy) from the OCDB
5049 //-----------------------------------------------------------------
5051 // check dead zones at z=0,+-7cm in the SPD
5052 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
5053 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5054 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
5055 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
5056 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5057 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
5058 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
5059 for (Int_t i=0; i<3; i++)
5060 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) return 1;
5063 // check bad zones from OCDB
5064 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
5066 if (idet<0) return 0;
5068 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
5070 // check if this detector is bad
5072 //printf("lay %d bad detector %d\n",ilayer,idet);
5077 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
5078 if (ilayer==0 || ilayer==1) { // ---------- SPD
5080 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
5082 detSizeFactorX *= 2.;
5083 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
5086 AliITSsegmentation *segm = (AliITSsegmentation*)fDetTypeRec->GetSegmentationModel(detType);
5087 if (detType==2) segm->SetLayer(ilayer+1);
5088 Float_t detSizeX = detSizeFactorX*segm->Dx();
5089 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
5091 // check if the road overlaps with bad chips
5093 LocalModuleCoord(ilayer,idet,track,xloc,zloc);
5094 Float_t zlocmin = zloc-dz;
5095 Float_t zlocmax = zloc+dz;
5096 Float_t xlocmin = xloc-dy;
5097 Float_t xlocmax = xloc+dy;
5098 Int_t chipsInRoad[100];
5100 if (TMath::Max(TMath::Abs(xlocmin),TMath::Abs(xlocmax))>0.5*detSizeX ||
5101 TMath::Max(TMath::Abs(zlocmin),TMath::Abs(zlocmax))>0.5*detSizeZ) return 0;
5102 //printf("lay %d det %d zmim zmax %f %f xmin xmax %f %f %f %f\n",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,segm->Dx(),segm->Dz());
5103 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
5104 //printf("lay %d nChipsInRoad %d\n",ilayer,nChipsInRoad);
5105 if (!nChipsInRoad) return 0;
5107 Bool_t anyBad=kFALSE,anyGood=kFALSE;
5108 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
5109 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
5110 //printf(" chip %d bad %d\n",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh]));
5111 if (det.IsChipBad(chipsInRoad[iCh])) {
5118 if (!anyGood) return 2; // all chips in road are bad
5120 if (anyBad) return 3; // at least a bad chip in road
5123 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
5124 || ilayer==4 || ilayer==5 // SSD
5125 || !noClusters) return 0;
5127 // There are no clusters in road: check if there is at least
5128 // a bad SPD pixel or SDD anode
5130 if(ilayer==1 || ilayer==3 || ilayer==5)
5131 idet += AliITSgeomTGeo::GetNLadders(ilayer)*AliITSgeomTGeo::GetNDetectors(ilayer);
5133 //if (fITSChannelStatus->AnyBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax)) return 3;
5135 if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
5139 //------------------------------------------------------------------------
5140 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
5141 AliITStrackMI *track,
5142 Float_t &xloc,Float_t &zloc) const {
5143 //-----------------------------------------------------------------
5144 // Gives position of track in local module ref. frame
5145 //-----------------------------------------------------------------
5150 if(idet<0) return kFALSE;
5152 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
5154 Int_t lad = Int_t(idet/ndet) + 1;
5156 Int_t det = idet - (lad-1)*ndet + 1;
5158 Double_t xyzGlob[3],xyzLoc[3];
5160 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
5161 // take into account the misalignment: xyz at real detector plane
5162 track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob);
5164 AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc);
5166 xloc = (Float_t)xyzLoc[0];
5167 zloc = (Float_t)xyzLoc[2];
5171 //------------------------------------------------------------------------
5172 Bool_t AliITStrackerMI::IsOKForPlaneEff(AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
5174 // Method to be optimized further:
5175 // Aim: decide whether a track can be used for PlaneEff evaluation
5176 // the decision is taken based on the track quality at the layer under study
5177 // no information on the clusters on this layer has to be used
5178 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
5179 // the cut is done on number of sigmas from the boundaries
5181 // Input: Actual track, layer [0,5] under study
5183 // Return: kTRUE if this is a good track
5185 // it will apply a pre-selection to obtain good quality tracks.
5186 // Here also you will have the possibility to put a control on the
5187 // impact point of the track on the basic block, in order to exclude border regions
5188 // this will be done by calling a proper method of the AliITSPlaneEff class.
5190 // input: AliITStrackMI* track, ilayer= layer number [0,5]
5191 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
5193 Int_t index[AliITSgeomTGeo::kNLayers];
5195 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
5197 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
5198 index[k]=clusters[k];
5202 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
5203 AliITSlayer &layer=fgLayers[ilayer];
5204 Double_t r=layer.GetR();
5205 AliITStrackMI tmp(*track);
5207 // require a minimal number of cluster in other layers and eventually clusters in closest layers
5209 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
5210 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
5211 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
5212 if (tmp.GetClIndex(lay)>0) ncl++;
5214 Bool_t nextout = kFALSE;
5215 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
5216 else nextout = ((tmp.GetClIndex(ilayer+1)>0)? kTRUE : kFALSE );
5217 Bool_t nextin = kFALSE;
5218 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
5219 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
5220 if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
5222 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
5223 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
5224 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
5225 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
5229 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
5230 Int_t idet=layer.FindDetectorIndex(phi,z);
5231 if(idet<0) { AliInfo(Form("cannot find detector"));
5234 // here check if it has good Chi Square.
5236 //propagate to the intersection with the detector plane
5237 const AliITSdetector &det=layer.GetDetector(idet);
5238 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
5242 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
5243 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5244 if(key>fPlaneEff->Nblock()) return kFALSE;
5245 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
5246 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
5248 // DEFINITION OF SEARCH ROAD FOR accepting a track
5250 //For the time being they are hard-wired, later on from AliITSRecoParam
5251 // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
5252 // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
5255 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5256 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5257 // done for RecPoints
5259 // exclude tracks at boundary between detectors
5260 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
5261 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
5262 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
5263 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
5264 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5266 if ( (locx-dx < blockXmn+boundaryWidth) ||
5267 (locx+dx > blockXmx-boundaryWidth) ||
5268 (locz-dz < blockZmn+boundaryWidth) ||
5269 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5272 //------------------------------------------------------------------------
5273 void AliITStrackerMI::UseTrackForPlaneEff(AliITStrackMI* track, Int_t ilayer) {
5275 // This Method has to be optimized! For the time-being it uses the same criteria
5276 // as those used in the search of extra clusters for overlapping modules.
5278 // Method Purpose: estabilish whether a track has produced a recpoint or not
5279 // in the layer under study (For Plane efficiency)
5281 // inputs: AliITStrackMI* track (pointer to a usable track)
5283 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5284 // traversed by this very track. In details:
5285 // - if a cluster can be associated to the track then call
5286 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5287 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5290 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5291 AliITSlayer &layer=fgLayers[ilayer];
5292 Double_t r=layer.GetR();
5293 AliITStrackMI tmp(*track);
5297 if (!tmp.GetPhiZat(r,phi,z)) return;
5298 Int_t idet=layer.FindDetectorIndex(phi,z);
5300 if(idet<0) { AliInfo(Form("cannot find detector"));
5304 //propagate to the intersection with the detector plane
5305 const AliITSdetector &det=layer.GetDetector(idet);
5306 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5310 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5312 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5313 TMath::Sqrt(tmp.GetSigmaZ2() +
5314 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5315 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5316 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5317 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5318 TMath::Sqrt(tmp.GetSigmaY2() +
5319 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5320 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5321 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5323 // road in global (rphi,z) [i.e. in tracking ref. system]
5324 Double_t zmin = tmp.GetZ() - dz;
5325 Double_t zmax = tmp.GetZ() + dz;
5326 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5327 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5329 // select clusters in road
5330 layer.SelectClusters(zmin,zmax,ymin,ymax);
5332 // Define criteria for track-cluster association
5333 Double_t msz = tmp.GetSigmaZ2() +
5334 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5335 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5336 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5337 Double_t msy = tmp.GetSigmaY2() +
5338 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5339 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5340 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5341 if (tmp.GetConstrain()) {
5342 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5343 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5345 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5346 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5348 msz = 1./msz; // 1/RoadZ^2
5349 msy = 1./msy; // 1/RoadY^2
5352 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5354 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5355 //Double_t tolerance=0.2;
5356 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5357 idetc = cl->GetDetectorIndex();
5358 if(idet!=idetc) continue;
5359 //Int_t ilay = cl->GetLayer();
5361 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5362 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5364 Double_t chi2=tmp.GetPredictedChi2(cl);
5365 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5369 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5371 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5372 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5373 if(key>fPlaneEff->Nblock()) return;
5374 Bool_t found=kFALSE;
5377 while ((cl=layer.GetNextCluster(clidx))!=0) {
5378 idetc = cl->GetDetectorIndex();
5379 if(idet!=idetc) continue;
5380 // here real control to see whether the cluster can be associated to the track.
5381 // cluster not associated to track
5382 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5383 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
5384 // calculate track-clusters chi2
5385 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5386 // in particular, the error associated to the cluster
5387 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5389 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5391 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5392 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5393 // track->SetExtraModule(ilayer,idetExtra);
5395 if(!fPlaneEff->UpDatePlaneEff(found,key))
5396 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5397 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5398 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
5399 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5400 Int_t cltype[2]={-999,-999};
5404 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
5405 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
5408 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5409 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5410 cltype[0]=layer.GetCluster(ci)->GetNy();
5411 cltype[1]=layer.GetCluster(ci)->GetNz();
5413 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5414 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
5415 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5416 // It is computed properly by calling the method
5417 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5419 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5420 //if (tmp.PropagateTo(x,0.,0.)) {
5421 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5422 AliCluster c(*layer.GetCluster(ci));
5423 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5424 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5425 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
5426 clu[2]=TMath::Sqrt(c.GetSigmaY2());
5427 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5430 fPlaneEff->FillHistos(key,found,tr,clu,cltype);