1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------------------------
18 // Implementation of the ITS tracker class
19 // It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
20 // and fills with them the ESD
21 // Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
22 // Current support and development:
23 // Andrea Dainese, andrea.dainese@lnl.infn.it
24 // dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
25 // Params moved to AliITSRecoParam by: Andrea Dainese, INFN
26 // Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
27 //-------------------------------------------------------------------------
31 #include <TDatabasePDG.h>
34 #include <TTreeStream.h>
38 #include "AliGeomManager.h"
39 #include "AliITSPlaneEff.h"
40 #include "AliTrackPointArray.h"
41 #include "AliESDEvent.h"
42 #include "AliESDtrack.h"
44 #include "AliITSChannelStatus.h"
45 #include "AliITSDetTypeRec.h"
46 #include "AliITSRecPoint.h"
47 #include "AliITSRecPointContainer.h"
48 #include "AliITSgeomTGeo.h"
49 #include "AliITSReconstructor.h"
50 #include "AliITSClusterParam.h"
51 #include "AliITSsegmentation.h"
52 #include "AliITSCalibration.h"
53 #include "AliITSPlaneEffSPD.h"
54 #include "AliITSPlaneEffSDD.h"
55 #include "AliITSPlaneEffSSD.h"
56 #include "AliITSV0Finder.h"
57 #include "AliITStrackerMI.h"
58 #include "AliMathBase.h"
60 ClassImp(AliITStrackerMI)
62 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
64 AliITStrackerMI::AliITStrackerMI():AliTracker(),
74 fLastLayerToTrackTo(0),
77 fTrackingPhase("Default"),
83 fxTimesRhoPipeTrks(0),
84 fxOverX0ShieldTrks(0),
85 fxTimesRhoShieldTrks(0),
87 fxTimesRhoLayerTrks(0),
94 for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
95 for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
96 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 !");
135 fOriginal.SetOwner();
139 for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
140 Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
141 Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
143 Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
144 AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
145 Double_t poff=TMath::ATan2(y,x);
148 AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
149 Double_t r=TMath::Sqrt(x*x + y*y);
151 AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
152 r += TMath::Sqrt(x*x + y*y);
153 AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
154 r += TMath::Sqrt(x*x + y*y);
155 AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
156 r += TMath::Sqrt(x*x + y*y);
159 new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
161 for (Int_t j=1; j<nlad+1; j++) {
162 for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
163 TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
164 const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
166 Double_t txyz[3]={0.};
167 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
168 m.LocalToMaster(txyz,xyz);
169 r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
170 Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
172 if (phi<0) phi+=TMath::TwoPi();
173 else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
175 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
176 new(&det) AliITSdetector(r,phi);
177 // compute the real radius (with misalignment)
178 TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
180 xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
181 mmisal.LocalToMaster(txyz,xyz);
182 Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
183 det.SetRmisal(rmisal);
185 } // end loop on detectors
186 } // end loop on ladders
187 fForceSkippingOfLayer[i] = 0;
188 } // end loop on layers
191 fI=AliITSgeomTGeo::GetNLayers();
194 fConstraint[0]=1; fConstraint[1]=0;
196 Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
197 AliITSReconstructor::GetRecoParam()->GetYVdef(),
198 AliITSReconstructor::GetRecoParam()->GetZVdef()};
199 Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
200 AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
201 AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
202 SetVertex(xyzVtx,ersVtx);
204 fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
205 for (Int_t i=0;i<100000;i++){
206 fBestTrackIndex[i]=0;
209 // store positions of centre of SPD modules (in z)
210 // The convetion is that fSPDdetzcentre is ordered from -z to +z
212 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
213 if (tr[2]<0) { // old geom (up to v5asymmPPR)
214 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
215 fSPDdetzcentre[0] = tr[2];
216 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
217 fSPDdetzcentre[1] = tr[2];
218 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
219 fSPDdetzcentre[2] = tr[2];
220 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
221 fSPDdetzcentre[3] = tr[2];
222 } else { // new geom (from v11Hybrid)
223 AliITSgeomTGeo::GetTranslation(1,1,4,tr);
224 fSPDdetzcentre[0] = tr[2];
225 AliITSgeomTGeo::GetTranslation(1,1,3,tr);
226 fSPDdetzcentre[1] = tr[2];
227 AliITSgeomTGeo::GetTranslation(1,1,2,tr);
228 fSPDdetzcentre[2] = tr[2];
229 AliITSgeomTGeo::GetTranslation(1,1,1,tr);
230 fSPDdetzcentre[3] = tr[2];
233 fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
234 if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
235 AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
239 for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
240 for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
242 if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
243 fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
245 // only for plane efficiency evaluation
246 if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
247 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
248 Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
249 if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
250 AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
251 if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
252 else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
253 else fPlaneEff = new AliITSPlaneEffSSD();
254 if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
255 if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
256 if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
259 //------------------------------------------------------------------------
260 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
262 fBestTrack(tracker.fBestTrack),
263 fTrackToFollow(tracker.fTrackToFollow),
264 fTrackHypothesys(tracker.fTrackHypothesys),
265 fBestHypothesys(tracker.fBestHypothesys),
266 fOriginal(tracker.fOriginal),
267 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
268 fPass(tracker.fPass),
269 fAfterV0(tracker.fAfterV0),
270 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
271 fCoefficients(tracker.fCoefficients),
273 fTrackingPhase(tracker.fTrackingPhase),
274 fUseTGeo(tracker.fUseTGeo),
275 fNtracks(tracker.fNtracks),
276 fxOverX0Pipe(tracker.fxOverX0Pipe),
277 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
279 fxTimesRhoPipeTrks(0),
280 fxOverX0ShieldTrks(0),
281 fxTimesRhoShieldTrks(0),
282 fxOverX0LayerTrks(0),
283 fxTimesRhoLayerTrks(0),
284 fDebugStreamer(tracker.fDebugStreamer),
285 fITSChannelStatus(tracker.fITSChannelStatus),
286 fkDetTypeRec(tracker.fkDetTypeRec),
287 fPlaneEff(tracker.fPlaneEff) {
289 fOriginal.SetOwner();
292 fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
295 fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
296 fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
299 fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
300 fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
303 //------------------------------------------------------------------------
304 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
305 //Assignment operator
306 this->~AliITStrackerMI();
307 new(this) AliITStrackerMI(tracker);
310 //------------------------------------------------------------------------
311 AliITStrackerMI::~AliITStrackerMI()
316 if (fCoefficients) delete [] fCoefficients;
317 DeleteTrksMaterialLUT();
318 if (fDebugStreamer) {
319 //fDebugStreamer->Close();
320 delete fDebugStreamer;
322 if(fITSChannelStatus) delete fITSChannelStatus;
323 if(fPlaneEff) delete fPlaneEff;
325 //------------------------------------------------------------------------
326 void AliITStrackerMI::ReadBadFromDetTypeRec() {
327 //--------------------------------------------------------------------
328 //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
330 //--------------------------------------------------------------------
332 if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
334 Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
336 if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
339 if(fITSChannelStatus) delete fITSChannelStatus;
340 fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
342 // ITS detectors and chips
343 Int_t i=0,j=0,k=0,ndet=0;
344 for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
345 Int_t nBadDetsPerLayer=0;
346 ndet=AliITSgeomTGeo::GetNDetectors(i);
347 for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
348 for (k=1; k<ndet+1; k++) {
349 AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
350 det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
351 if(det.IsBad()) {nBadDetsPerLayer++;}
352 } // end loop on detectors
353 } // end loop on ladders
354 AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
355 } // end loop on layers
359 //------------------------------------------------------------------------
360 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
361 //--------------------------------------------------------------------
362 //This function loads ITS clusters
363 //--------------------------------------------------------------------
365 TClonesArray *clusters = NULL;
366 AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
367 clusters=rpcont->FetchClusters(0,cTree);
368 if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
369 AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
372 Int_t i=0,j=0,ndet=0;
374 for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
375 ndet=fgLayers[i].GetNdetectors();
376 Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
377 for (; j<jmax; j++) {
378 // if (!cTree->GetEvent(j)) continue;
379 clusters = rpcont->UncheckedGetClusters(j);
380 if(!clusters)continue;
381 Int_t ncl=clusters->GetEntriesFast();
382 SignDeltas(clusters,GetZ());
385 AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
386 detector=c->GetDetectorIndex();
388 if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
390 Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
392 AliWarning(Form("Too many clusters on layer %d!",i));
397 // add dead zone "virtual" cluster in SPD, if there is a cluster within
398 // zwindow cm from the dead zone
399 // This method assumes that fSPDdetzcentre is ordered from -z to +z
400 if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
401 for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
402 Int_t lab[4] = {0,0,0,detector};
403 Int_t info[3] = {0,0,i};
404 Float_t q = 0.; // this identifies virtual clusters
405 Float_t hit[5] = {xdead,
407 AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
408 AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
410 Bool_t local = kTRUE;
411 Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
412 hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
413 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
414 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
415 hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
416 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
417 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
418 hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
419 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
420 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
421 hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
422 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
423 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
424 hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
425 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
426 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
427 hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
428 if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
429 fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
431 } // "virtual" clusters in SPD
435 fgLayers[i].ResetRoad(); //road defined by the cluster density
436 fgLayers[i].SortClusters();
439 // check whether we have to skip some layers
440 SetForceSkippingOfLayer();
444 //------------------------------------------------------------------------
445 void AliITStrackerMI::UnloadClusters() {
446 //--------------------------------------------------------------------
447 //This function unloads ITS clusters
448 //--------------------------------------------------------------------
449 for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
451 //------------------------------------------------------------------------
452 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
453 //--------------------------------------------------------------------
454 // Publishes all pointers to clusters known to the tracker into the
455 // passed object array.
456 // The ownership is not transfered - the caller is not expected to delete
458 //--------------------------------------------------------------------
460 for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
461 for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
462 AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
469 //------------------------------------------------------------------------
470 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
471 //--------------------------------------------------------------------
472 // Correction for the material between the TPC and the ITS
473 //--------------------------------------------------------------------
474 if (t->GetX() > AliITSRecoParam::Getriw()) { // inward direction
475 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
476 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
477 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
478 } else if (t->GetX() < AliITSRecoParam::Getrs()) { // outward direction
479 if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1)) return 0;// ITS screen
480 if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
481 if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
483 printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
489 //------------------------------------------------------------------------
490 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
491 //--------------------------------------------------------------------
492 // This functions reconstructs ITS tracks
493 // The clusters must be already loaded !
494 //--------------------------------------------------------------------
496 AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
498 fTrackingPhase="Clusters2Tracks";
500 TObjArray itsTracks(15000);
502 fEsd = event; // store pointer to the esd
504 // temporary (for cosmics)
505 if(event->GetVertex()) {
506 TString title = event->GetVertex()->GetTitle();
507 if(title.Contains("cosmics")) {
508 Double_t xyz[3]={GetX(),GetY(),GetZ()};
509 Double_t exyz[3]={0.1,0.1,0.1};
515 {/* Read ESD tracks */
516 Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
517 Int_t nentr=event->GetNumberOfTracks();
519 // Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
521 AliESDtrack *esd=event->GetTrack(nentr);
522 // ---- for debugging:
523 //if(TMath::Abs(esd->GetX()-83.65)<0.1) { FILE *f=fopen("tpc.dat","a"); fprintf(f,"%f %f %f %f %f %f\n",(Float_t)event->GetEventNumberInFile(),(Float_t)TMath::Abs(esd->GetLabel()),(Float_t)esd->GetX(),(Float_t)esd->GetY(),(Float_t)esd->GetZ(),(Float_t)esd->Pt()); fclose(f); }
525 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
526 if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
527 if (esd->GetStatus()&AliESDtrack::kITSin) continue;
528 if (esd->GetKinkIndex(0)>0) continue; //kink daughter
529 AliITStrackMI *t = new AliITStrackMI(*esd);
530 t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP()); //I.B.
531 Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
534 // look at the ESD mass hypothesys !
535 if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
537 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
539 if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
540 //track - can be V0 according to TPC
542 if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
546 if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
550 if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
554 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
559 t->SetReconstructed(kFALSE);
560 itsTracks.AddLast(t);
561 fOriginal.AddLast(t);
563 } /* End Read ESD tracks */
567 Int_t nentr=itsTracks.GetEntriesFast();
568 fTrackHypothesys.Expand(nentr);
569 fBestHypothesys.Expand(nentr);
570 MakeCoefficients(nentr);
571 if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
573 // THE TWO TRACKING PASSES
574 for (fPass=0; fPass<2; fPass++) {
575 Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
576 for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
577 AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
578 if (t==0) continue; //this track has been already tracked
579 //cout<<"========== "<<fPass<<" "<<fCurrentEsdTrack<<" =========\n";
580 if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue; //this track was already "succesfully" reconstructed
581 Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz); //I.B.
582 if (fConstraint[fPass]) {
583 if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
584 TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
587 Int_t tpcLabel=t->GetLabel(); //save the TPC track label
588 AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
590 ResetTrackToFollow(*t);
593 FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
596 SortTrackHypothesys(fCurrentEsdTrack,20,0); //MI change
598 AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
599 if (!besttrack) continue;
600 besttrack->SetLabel(tpcLabel);
601 // besttrack->CookdEdx();
603 besttrack->SetFakeRatio(1.);
604 CookLabel(besttrack,0.); //For comparison only
605 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
607 if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue; //to be tracked also without vertex constrain
609 t->SetReconstructed(kTRUE);
611 AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
613 GetBestHypothesysMIP(itsTracks);
614 } // end loop on the two tracking passes
616 if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
617 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
622 Int_t entries = fTrackHypothesys.GetEntriesFast();
623 for (Int_t ientry=0; ientry<entries; ientry++) {
624 TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
625 if (array) array->Delete();
626 delete fTrackHypothesys.RemoveAt(ientry);
629 fTrackHypothesys.Delete();
630 entries = fBestHypothesys.GetEntriesFast();
631 for (Int_t ientry=0; ientry<entries; ientry++) {
632 TObjArray * array =(TObjArray*)fBestHypothesys.UncheckedAt(ientry);
633 if (array) array->Delete();
634 delete fBestHypothesys.RemoveAt(ientry);
636 fBestHypothesys.Delete();
638 delete [] fCoefficients;
640 DeleteTrksMaterialLUT();
642 AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
644 fTrackingPhase="Default";
648 //------------------------------------------------------------------------
649 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
650 //--------------------------------------------------------------------
651 // This functions propagates reconstructed ITS tracks back
652 // The clusters must be loaded !
653 //--------------------------------------------------------------------
654 fTrackingPhase="PropagateBack";
655 Int_t nentr=event->GetNumberOfTracks();
656 // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
659 for (Int_t i=0; i<nentr; i++) {
660 AliESDtrack *esd=event->GetTrack(i);
662 if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
663 if (esd->GetStatus()&AliESDtrack::kITSout) continue;
665 AliITStrackMI *t = new AliITStrackMI(*esd);
667 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
669 ResetTrackToFollow(*t);
672 // propagate to vertex [SR, GSI 17.02.2003]
673 // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
674 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
675 if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
676 fTrackToFollow.StartTimeIntegral();
677 // from vertex to outside pipe
678 CorrectForPipeMaterial(&fTrackToFollow,"outward");
680 // Start time integral and add distance from current position to vertex
681 Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
682 fTrackToFollow.GetXYZ(xyzTrk);
684 for (Int_t icoord=0; icoord<3; icoord++) {
685 Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
688 fTrackToFollow.StartTimeIntegral();
689 fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
691 fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
692 if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
693 if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
697 fTrackToFollow.SetLabel(t->GetLabel());
698 //fTrackToFollow.CookdEdx();
699 CookLabel(&fTrackToFollow,0.); //For comparison only
700 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
701 //UseClusters(&fTrackToFollow);
707 AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
709 fTrackingPhase="Default";
713 //------------------------------------------------------------------------
714 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
715 //--------------------------------------------------------------------
716 // This functions refits ITS tracks using the
717 // "inward propagated" TPC tracks
718 // The clusters must be loaded !
719 //--------------------------------------------------------------------
720 fTrackingPhase="RefitInward";
722 if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
724 Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
725 if(!doExtra) AliDebug(2,"Do not search for extra clusters");
727 Int_t nentr=event->GetNumberOfTracks();
728 // Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
731 for (Int_t i=0; i<nentr; i++) {
732 AliESDtrack *esd=event->GetTrack(i);
734 if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
735 if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
736 if (esd->GetStatus()&AliESDtrack::kTPCout)
737 if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
739 AliITStrackMI *t = new AliITStrackMI(*esd);
741 t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
742 if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
747 ResetTrackToFollow(*t);
748 fTrackToFollow.ResetClusters();
750 // ITS standalone tracks
751 if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
752 fTrackToFollow.ResetCovariance(10.);
753 // protection for loopers that can have parameters screwed up
754 if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
755 TMath::Abs(fTrackToFollow.GetZ())>1000.) {
762 Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
763 AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
765 AliDebug(2,Form("Refit LABEL %d %d",t->GetLabel(),t->GetNumberOfClusters()));
766 if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
767 AliDebug(2," refit OK");
768 fTrackToFollow.SetLabel(t->GetLabel());
769 // fTrackToFollow.CookdEdx();
770 CookdEdx(&fTrackToFollow);
772 CookLabel(&fTrackToFollow,0.0); //For comparison only
775 if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
776 fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
777 AliESDtrack *esdTrack =fTrackToFollow.GetESDtrack();
778 //printf(" %d\n",esdTrack->GetITSModuleIndex(0));
779 //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
780 Double_t r[3]={0.,0.,0.};
782 esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
789 AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
791 fTrackingPhase="Default";
795 //------------------------------------------------------------------------
796 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
797 //--------------------------------------------------------------------
798 // Return pointer to a given cluster
799 //--------------------------------------------------------------------
800 Int_t l=(index & 0xf0000000) >> 28;
801 Int_t c=(index & 0x0fffffff) >> 00;
802 return fgLayers[l].GetCluster(c);
804 //------------------------------------------------------------------------
805 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
806 //--------------------------------------------------------------------
807 // Get track space point with index i
808 //--------------------------------------------------------------------
810 Int_t l=(index & 0xf0000000) >> 28;
811 Int_t c=(index & 0x0fffffff) >> 00;
812 AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
813 Int_t idet = cl->GetDetectorIndex();
817 cl->GetGlobalXYZ(xyz);
818 cl->GetGlobalCov(cov);
820 p.SetCharge(cl->GetQ());
821 p.SetDriftTime(cl->GetDriftTime());
822 p.SetChargeRatio(cl->GetChargeRatio());
823 p.SetClusterType(cl->GetClusterType());
824 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
827 iLayer = AliGeomManager::kSPD1;
830 iLayer = AliGeomManager::kSPD2;
833 iLayer = AliGeomManager::kSDD1;
836 iLayer = AliGeomManager::kSDD2;
839 iLayer = AliGeomManager::kSSD1;
842 iLayer = AliGeomManager::kSSD2;
845 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
848 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
849 p.SetVolumeID((UShort_t)volid);
852 //------------------------------------------------------------------------
853 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
854 AliTrackPoint& p, const AliESDtrack *t) {
855 //--------------------------------------------------------------------
856 // Get track space point with index i
857 // (assign error estimated during the tracking)
858 //--------------------------------------------------------------------
860 Int_t l=(index & 0xf0000000) >> 28;
861 Int_t c=(index & 0x0fffffff) >> 00;
862 const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
863 Int_t idet = cl->GetDetectorIndex();
865 const AliITSdetector &det=fgLayers[l].GetDetector(idet);
867 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
869 detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
870 detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
871 Double_t alpha = t->GetAlpha();
872 Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
873 Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
874 phi += alpha-det.GetPhi();
875 Float_t tgphi = TMath::Tan(phi);
877 Float_t tgl = t->GetTgl(); // tgl about const along track
878 Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
880 Float_t errtrky,errtrkz,covyz;
881 Bool_t addMisalErr=kFALSE;
882 AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
886 cl->GetGlobalXYZ(xyz);
887 // cl->GetGlobalCov(cov);
888 Float_t pos[3] = {0.,0.,0.};
889 AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
890 tmpcl.GetGlobalCov(cov);
893 p.SetCharge(cl->GetQ());
894 p.SetDriftTime(cl->GetDriftTime());
895 p.SetChargeRatio(cl->GetChargeRatio());
896 p.SetClusterType(cl->GetClusterType());
898 AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
901 iLayer = AliGeomManager::kSPD1;
904 iLayer = AliGeomManager::kSPD2;
907 iLayer = AliGeomManager::kSDD1;
910 iLayer = AliGeomManager::kSDD2;
913 iLayer = AliGeomManager::kSSD1;
916 iLayer = AliGeomManager::kSSD2;
919 AliWarning(Form("Wrong layer index in ITS (%d) !",l));
922 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
924 p.SetVolumeID((UShort_t)volid);
927 //------------------------------------------------------------------------
928 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
930 //--------------------------------------------------------------------
931 // Follow prolongation tree
932 //--------------------------------------------------------------------
934 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
935 Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
938 AliESDtrack * esd = otrack->GetESDtrack();
939 if (esd->GetV0Index(0)>0) {
940 // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
941 // mapping of ESD track is different as ITS track in Containers
942 // Need something more stable
943 // Indexes are set back again to the ESD track indexes in UpdateTPCV0
944 for (Int_t i=0;i<3;i++){
945 Int_t index = esd->GetV0Index(i);
947 AliESDv0 * vertex = fEsd->GetV0(index);
948 if (vertex->GetStatus()<0) continue; // rejected V0
950 if (esd->GetSign()>0) {
951 vertex->SetIndex(0,esdindex);
953 vertex->SetIndex(1,esdindex);
957 TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
959 bestarray = new TObjArray(5);
960 bestarray->SetOwner();
961 fBestHypothesys.AddAt(bestarray,esdindex);
965 //setup tree of the prolongations
967 static AliITStrackMI tracks[7][100];
968 AliITStrackMI *currenttrack;
969 static AliITStrackMI currenttrack1;
970 static AliITStrackMI currenttrack2;
971 static AliITStrackMI backuptrack;
973 Int_t nindexes[7][100];
974 Float_t normalizedchi2[100];
975 for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
976 otrack->SetNSkipped(0);
977 new (&(tracks[6][0])) AliITStrackMI(*otrack);
979 for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
980 Int_t modstatus = 1; // found
984 // follow prolongations
985 for (Int_t ilayer=5; ilayer>=0; ilayer--) {
986 AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
989 AliITSlayer &layer=fgLayers[ilayer];
990 Double_t r = layer.GetR();
996 for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
998 if (ntracks[ilayer]>=100) break;
999 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1000 if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1001 if (ntracks[ilayer]>15+ilayer){
1002 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1003 if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1006 new(¤ttrack1) AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1008 // material between SSD and SDD, SDD and SPD
1010 if(!CorrectForShieldMaterial(¤ttrack1,"SDD","inward")) continue;
1012 if(!CorrectForShieldMaterial(¤ttrack1,"SPD","inward")) continue;
1016 if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1017 Int_t idet=layer.FindDetectorIndex(phi,z);
1019 Double_t trackGlobXYZ1[3];
1020 if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1022 // Get the budget to the primary vertex for the current track being prolonged
1023 Double_t budgetToPrimVertex = GetEffectiveThickness();
1025 // check if we allow a prolongation without point
1026 Int_t skip = CheckSkipLayer(¤ttrack1,ilayer,idet);
1028 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1029 // propagate to the layer radius
1030 Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1031 if(!vtrack->Propagate(xToGo)) continue;
1032 // apply correction for material of the current layer
1033 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1034 vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1035 vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1036 vtrack->SetClIndex(ilayer,-1);
1037 modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1038 if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1039 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1041 if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1046 // track outside layer acceptance in z
1047 if (idet<0) continue;
1049 //propagate to the intersection with the detector plane
1050 const AliITSdetector &det=layer.GetDetector(idet);
1051 new(¤ttrack2) AliITStrackMI(currenttrack1);
1052 if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1053 if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1054 currenttrack1.SetDetectorIndex(idet);
1055 currenttrack2.SetDetectorIndex(idet);
1056 if(!LocalModuleCoord(ilayer,idet,¤ttrack1,xloc,zloc)) continue; // local module coords
1059 // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1061 // road in global (rphi,z) [i.e. in tracking ref. system]
1062 Double_t zmin,zmax,ymin,ymax;
1063 if (!ComputeRoad(¤ttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1065 // select clusters in road
1066 layer.SelectClusters(zmin,zmax,ymin,ymax);
1067 //********************
1069 // Define criteria for track-cluster association
1070 Double_t msz = currenttrack1.GetSigmaZ2() +
1071 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1072 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1073 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1074 Double_t msy = currenttrack1.GetSigmaY2() +
1075 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1076 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1077 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1079 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1080 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1082 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1083 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1085 msz = 1./msz; // 1/RoadZ^2
1086 msy = 1./msy; // 1/RoadY^2
1090 // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1092 const AliITSRecPoint *cl=0;
1094 Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1095 Bool_t deadzoneSPD=kFALSE;
1096 currenttrack = ¤ttrack1;
1098 // check if the road contains a dead zone
1099 Bool_t noClusters = kFALSE;
1100 if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1101 if (noClusters) AliDebug(2,"no clusters in road");
1102 Double_t dz=0.5*(zmax-zmin);
1103 Double_t dy=0.5*(ymax-ymin);
1104 Int_t dead = CheckDeadZone(¤ttrack1,ilayer,idet,dz,dy,noClusters);
1105 if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1106 // create a prolongation without clusters (check also if there are no clusters in the road)
1109 AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1110 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1111 updatetrack->SetClIndex(ilayer,-1);
1113 modstatus = 5; // no cls in road
1114 } else if (dead==1) {
1115 modstatus = 7; // holes in z in SPD
1116 } else if (dead==2 || dead==3 || dead==4) {
1117 modstatus = 2; // dead from OCDB
1119 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1120 // apply correction for material of the current layer
1121 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1122 if (constrain) { // apply vertex constrain
1123 updatetrack->SetConstrain(constrain);
1124 Bool_t isPrim = kTRUE;
1125 if (ilayer<4) { // check that it's close to the vertex
1126 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1127 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1128 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1129 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1130 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1132 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1134 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1136 if (dead==1) { // dead zone at z=0,+-7cm in SPD
1137 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1139 } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1140 updatetrack->SetDeadZoneProbability(ilayer,1.);
1141 } else if (dead==4) { // at least a single dead channel from OCDB
1142 updatetrack->SetDeadZoneProbability(ilayer,0.);
1149 // loop over clusters in the road
1150 while ((cl=layer.GetNextCluster(clidx))!=0) {
1151 if (ntracks[ilayer]>95) break; //space for skipped clusters
1152 Bool_t changedet =kFALSE;
1153 if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1154 Int_t idetc=cl->GetDetectorIndex();
1156 if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1157 // take into account misalignment (bring track to real detector plane)
1158 Double_t xTrOrig = currenttrack->GetX();
1159 if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1160 // a first cut on track-cluster distance
1161 if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1162 (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1163 { // cluster not associated to track
1164 AliDebug(2,"not associated");
1165 // MvL: added here as well
1166 // bring track back to ideal detector plane
1167 currenttrack->Propagate(xTrOrig);
1170 // bring track back to ideal detector plane
1171 if (!currenttrack->Propagate(xTrOrig)) continue;
1172 } else { // have to move track to cluster's detector
1173 const AliITSdetector &detc=layer.GetDetector(idetc);
1174 // a first cut on track-cluster distance
1176 if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1177 if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1178 (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1179 continue; // cluster not associated to track
1181 new (&backuptrack) AliITStrackMI(currenttrack2);
1183 currenttrack =¤ttrack2;
1184 if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1185 new (currenttrack) AliITStrackMI(backuptrack);
1189 currenttrack->SetDetectorIndex(idetc);
1190 // Get again the budget to the primary vertex
1191 // for the current track being prolonged, if had to change detector
1192 //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1195 // calculate track-clusters chi2
1196 chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1198 AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1199 if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1200 if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1201 if (ntracks[ilayer]>=100) continue;
1202 AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1203 updatetrack->SetClIndex(ilayer,-1);
1204 if (changedet) new (¤ttrack2) AliITStrackMI(backuptrack);
1206 if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1207 if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1208 AliDebug(2,"update failed");
1211 updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1212 modstatus = 1; // found
1213 } else { // virtual cluster in dead zone
1214 updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1215 updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1216 modstatus = 7; // holes in z in SPD
1220 Float_t xlocnewdet,zlocnewdet;
1221 if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1222 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1225 updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1227 if (cl->IsUsed()) updatetrack->IncrementNUsed();
1229 // apply correction for material of the current layer
1230 CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1232 if (constrain) { // apply vertex constrain
1233 updatetrack->SetConstrain(constrain);
1234 Bool_t isPrim = kTRUE;
1235 if (ilayer<4) { // check that it's close to the vertex
1236 updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1237 if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1238 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1239 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1240 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1242 if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1243 } //apply vertex constrain
1245 } // create new hypothesis
1247 AliDebug(2,"chi2 too large");
1250 } // loop over possible prolongations
1252 // allow one prolongation without clusters
1253 if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1254 AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1255 // apply correction for material of the current layer
1256 CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1257 vtrack->SetClIndex(ilayer,-1);
1258 modstatus = 3; // skipped
1259 vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1260 if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1261 vtrack->IncrementNSkipped();
1266 } // loop over tracks in layer ilayer+1
1268 //loop over track candidates for the current layer
1274 for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1275 normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1276 if (normalizedchi2[itrack] <
1277 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1281 if (constrain) { // constrain
1282 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1284 } else { // !constrain
1285 if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1290 // sort tracks by increasing normalized chi2
1291 TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1292 ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1293 if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1294 if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1295 } // end loop over layers
1299 // Now select tracks to be kept
1301 Int_t max = constrain ? 20 : 5;
1303 // tracks that reach layer 0 (SPD inner)
1304 for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1305 AliITStrackMI & track= tracks[0][nindexes[0][i]];
1306 if (track.GetNumberOfClusters()<2) continue;
1307 if (!constrain && track.GetNormChi2(0) >
1308 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1311 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1314 // tracks that reach layer 1 (SPD outer)
1315 for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1316 AliITStrackMI & track= tracks[1][nindexes[1][i]];
1317 if (track.GetNumberOfClusters()<4) continue;
1318 if (!constrain && track.GetNormChi2(1) >
1319 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1320 if (constrain) track.IncrementNSkipped();
1322 track.SetD(0,track.GetD(GetX(),GetY()));
1323 track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1324 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1325 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1328 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1331 // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1333 for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1334 AliITStrackMI & track= tracks[2][nindexes[2][i]];
1335 if (track.GetNumberOfClusters()<3) continue;
1336 if (!constrain && track.GetNormChi2(2) >
1337 AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1338 if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1340 track.SetD(0,track.GetD(GetX(),GetY()));
1341 track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1342 if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1343 track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1346 AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1352 // register best track of each layer - important for V0 finder
1354 for (Int_t ilayer=0;ilayer<5;ilayer++){
1355 if (ntracks[ilayer]==0) continue;
1356 AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1357 if (track.GetNumberOfClusters()<1) continue;
1358 CookLabel(&track,0);
1359 bestarray->AddAt(new AliITStrackMI(track),ilayer);
1363 // update TPC V0 information
1365 if (otrack->GetESDtrack()->GetV0Index(0)>0){
1366 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1367 for (Int_t i=0;i<3;i++){
1368 Int_t index = otrack->GetESDtrack()->GetV0Index(i);
1369 if (index==0) break;
1370 AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1371 if (vertex->GetStatus()<0) continue; // rejected V0
1373 if (otrack->GetSign()>0) {
1374 vertex->SetIndex(0,esdindex);
1377 vertex->SetIndex(1,esdindex);
1379 //find nearest layer with track info
1380 Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]); //I.B.
1381 Int_t nearestold = GetNearestLayer(xrp); //I.B.
1382 Int_t nearest = nearestold;
1383 for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1384 if (ntracks[nearest]==0){
1389 AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1390 if (nearestold<5&&nearest<5){
1391 Bool_t accept = track.GetNormChi2(nearest)<10;
1393 if (track.GetSign()>0) {
1394 vertex->SetParamP(track);
1395 vertex->Update(fprimvertex);
1396 //vertex->SetIndex(0,track.fESDtrack->GetID());
1397 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1399 vertex->SetParamN(track);
1400 vertex->Update(fprimvertex);
1401 //vertex->SetIndex(1,track.fESDtrack->GetID());
1402 if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1404 vertex->SetStatus(vertex->GetStatus()+1);
1406 //vertex->SetStatus(-2); // reject V0 - not enough clusters
1413 //------------------------------------------------------------------------
1414 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1416 //--------------------------------------------------------------------
1419 return fgLayers[layer];
1421 //------------------------------------------------------------------------
1422 AliITStrackerMI::AliITSlayer::AliITSlayer():
1452 //--------------------------------------------------------------------
1453 //default AliITSlayer constructor
1454 //--------------------------------------------------------------------
1455 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1456 fClusterWeight[i]=0;
1457 fClusterTracks[0][i]=-1;
1458 fClusterTracks[1][i]=-1;
1459 fClusterTracks[2][i]=-1;
1460 fClusterTracks[3][i]=-1;
1463 //------------------------------------------------------------------------
1464 AliITStrackerMI::AliITSlayer::
1465 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1494 //--------------------------------------------------------------------
1495 //main AliITSlayer constructor
1496 //--------------------------------------------------------------------
1497 fDetectors=new AliITSdetector[fNladders*fNdetectors];
1498 fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1500 //------------------------------------------------------------------------
1501 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1503 fPhiOffset(layer.fPhiOffset),
1504 fNladders(layer.fNladders),
1505 fZOffset(layer.fZOffset),
1506 fNdetectors(layer.fNdetectors),
1507 fDetectors(layer.fDetectors),
1512 fClustersCs(layer.fClustersCs),
1513 fClusterIndexCs(layer.fClusterIndexCs),
1517 fCurrentSlice(layer.fCurrentSlice),
1525 fAccepted(layer.fAccepted),
1527 fMaxSigmaClY(layer.fMaxSigmaClY),
1528 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1529 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1533 //------------------------------------------------------------------------
1534 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1535 //--------------------------------------------------------------------
1536 // AliITSlayer destructor
1537 //--------------------------------------------------------------------
1538 delete [] fDetectors;
1539 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1540 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1541 fClusterWeight[i]=0;
1542 fClusterTracks[0][i]=-1;
1543 fClusterTracks[1][i]=-1;
1544 fClusterTracks[2][i]=-1;
1545 fClusterTracks[3][i]=-1;
1548 //------------------------------------------------------------------------
1549 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1550 //--------------------------------------------------------------------
1551 // This function removes loaded clusters
1552 //--------------------------------------------------------------------
1553 for (Int_t i=0; i<fN; i++) delete fClusters[i];
1554 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1555 fClusterWeight[i]=0;
1556 fClusterTracks[0][i]=-1;
1557 fClusterTracks[1][i]=-1;
1558 fClusterTracks[2][i]=-1;
1559 fClusterTracks[3][i]=-1;
1565 //------------------------------------------------------------------------
1566 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1567 //--------------------------------------------------------------------
1568 // This function reset weights of the clusters
1569 //--------------------------------------------------------------------
1570 for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1571 fClusterWeight[i]=0;
1572 fClusterTracks[0][i]=-1;
1573 fClusterTracks[1][i]=-1;
1574 fClusterTracks[2][i]=-1;
1575 fClusterTracks[3][i]=-1;
1577 for (Int_t i=0; i<fN;i++) {
1578 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1579 if (cl&&cl->IsUsed()) cl->Use();
1583 //------------------------------------------------------------------------
1584 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1585 //--------------------------------------------------------------------
1586 // This function calculates the road defined by the cluster density
1587 //--------------------------------------------------------------------
1589 for (Int_t i=0; i<fN; i++) {
1590 if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1592 if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1594 //------------------------------------------------------------------------
1595 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1596 //--------------------------------------------------------------------
1597 //This function adds a cluster to this layer
1598 //--------------------------------------------------------------------
1599 if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1605 AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1607 Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1608 Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1609 if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1610 if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1611 if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1612 if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1615 if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1616 if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1617 if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1618 if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1622 //------------------------------------------------------------------------
1623 void AliITStrackerMI::AliITSlayer::SortClusters()
1628 AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1629 Float_t *z = new Float_t[fN];
1630 Int_t * index = new Int_t[fN];
1632 fMaxSigmaClY=0.; //AD
1633 fMaxSigmaClZ=0.; //AD
1635 for (Int_t i=0;i<fN;i++){
1636 z[i] = fClusters[i]->GetZ();
1637 // save largest errors in y and z for this layer
1638 fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1639 fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1641 TMath::Sort(fN,z,index,kFALSE);
1642 for (Int_t i=0;i<fN;i++){
1643 clusters[i] = fClusters[index[i]];
1646 for (Int_t i=0;i<fN;i++){
1647 fClusters[i] = clusters[i];
1648 fZ[i] = fClusters[i]->GetZ();
1649 AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1650 Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1651 if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1661 for (Int_t i=0;i<fN;i++){
1662 if (fY[i]<fYB[0]) fYB[0]=fY[i];
1663 if (fY[i]>fYB[1]) fYB[1]=fY[i];
1664 fClusterIndex[i] = i;
1668 fDy5 = (fYB[1]-fYB[0])/5.;
1669 fDy10 = (fYB[1]-fYB[0])/10.;
1670 fDy20 = (fYB[1]-fYB[0])/20.;
1671 for (Int_t i=0;i<6;i++) fN5[i] =0;
1672 for (Int_t i=0;i<11;i++) fN10[i]=0;
1673 for (Int_t i=0;i<21;i++) fN20[i]=0;
1675 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;}
1676 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;}
1677 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;}
1680 for (Int_t i=0;i<fN;i++)
1681 for (Int_t irot=-1;irot<=1;irot++){
1682 Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1684 for (Int_t slice=0; slice<6;slice++){
1685 if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1686 fClusters5[slice][fN5[slice]] = fClusters[i];
1687 fY5[slice][fN5[slice]] = curY;
1688 fZ5[slice][fN5[slice]] = fZ[i];
1689 fClusterIndex5[slice][fN5[slice]]=i;
1694 for (Int_t slice=0; slice<11;slice++){
1695 if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1696 fClusters10[slice][fN10[slice]] = fClusters[i];
1697 fY10[slice][fN10[slice]] = curY;
1698 fZ10[slice][fN10[slice]] = fZ[i];
1699 fClusterIndex10[slice][fN10[slice]]=i;
1704 for (Int_t slice=0; slice<21;slice++){
1705 if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1706 fClusters20[slice][fN20[slice]] = fClusters[i];
1707 fY20[slice][fN20[slice]] = curY;
1708 fZ20[slice][fN20[slice]] = fZ[i];
1709 fClusterIndex20[slice][fN20[slice]]=i;
1716 // consistency check
1718 for (Int_t i=0;i<fN-1;i++){
1724 for (Int_t slice=0;slice<21;slice++)
1725 for (Int_t i=0;i<fN20[slice]-1;i++){
1726 if (fZ20[slice][i]>fZ20[slice][i+1]){
1733 //------------------------------------------------------------------------
1734 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1735 //--------------------------------------------------------------------
1736 // This function returns the index of the nearest cluster
1737 //--------------------------------------------------------------------
1740 if (fCurrentSlice<0) {
1749 if (ncl==0) return 0;
1750 Int_t b=0, e=ncl-1, m=(b+e)/2;
1751 for (; b<e; m=(b+e)/2) {
1752 // if (z > fClusters[m]->GetZ()) b=m+1;
1753 if (z > zcl[m]) b=m+1;
1758 //------------------------------------------------------------------------
1759 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 {
1760 //--------------------------------------------------------------------
1761 // This function computes the rectangular road for this track
1762 //--------------------------------------------------------------------
1765 AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1766 // take into account the misalignment: propagate track to misaligned detector plane
1767 if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1769 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1770 TMath::Sqrt(track->GetSigmaZ2() +
1771 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1772 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1773 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1774 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1775 TMath::Sqrt(track->GetSigmaY2() +
1776 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1777 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1778 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1780 // track at boundary between detectors, enlarge road
1781 Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1782 if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1783 (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1784 (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1785 (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1786 Float_t tgl = TMath::Abs(track->GetTgl());
1787 if (tgl > 1.) tgl=1.;
1788 Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1789 dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1790 Float_t snp = TMath::Abs(track->GetSnp());
1791 if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1792 dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1795 // add to the road a term (up to 2-3 mm) to deal with misalignments
1796 dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1797 dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1799 Double_t r = fgLayers[ilayer].GetR();
1800 zmin = track->GetZ() - dz;
1801 zmax = track->GetZ() + dz;
1802 ymin = track->GetY() + r*det.GetPhi() - dy;
1803 ymax = track->GetY() + r*det.GetPhi() + dy;
1805 // bring track back to idead detector plane
1806 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1810 //------------------------------------------------------------------------
1811 void AliITStrackerMI::AliITSlayer::
1812 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1813 //--------------------------------------------------------------------
1814 // This function sets the "window"
1815 //--------------------------------------------------------------------
1817 Double_t circle=2*TMath::Pi()*fR;
1823 // enlarge road in y by maximum cluster error on this layer (3 sigma)
1824 fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1825 fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1826 fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1827 fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1829 Float_t ymiddle = (fYmin+fYmax)*0.5;
1830 if (ymiddle<fYB[0]) {
1831 fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1832 } else if (ymiddle>fYB[1]) {
1833 fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1839 fClustersCs = fClusters;
1840 fClusterIndexCs = fClusterIndex;
1846 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1847 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1848 if (slice<0) slice=0;
1849 if (slice>20) slice=20;
1850 Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1852 fCurrentSlice=slice;
1853 fClustersCs = fClusters20[fCurrentSlice];
1854 fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1855 fYcs = fY20[fCurrentSlice];
1856 fZcs = fZ20[fCurrentSlice];
1857 fNcs = fN20[fCurrentSlice];
1862 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1863 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1864 if (slice<0) slice=0;
1865 if (slice>10) slice=10;
1866 Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1868 fCurrentSlice=slice;
1869 fClustersCs = fClusters10[fCurrentSlice];
1870 fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1871 fYcs = fY10[fCurrentSlice];
1872 fZcs = fZ10[fCurrentSlice];
1873 fNcs = fN10[fCurrentSlice];
1878 if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1879 Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1880 if (slice<0) slice=0;
1881 if (slice>5) slice=5;
1882 Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1884 fCurrentSlice=slice;
1885 fClustersCs = fClusters5[fCurrentSlice];
1886 fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1887 fYcs = fY5[fCurrentSlice];
1888 fZcs = fZ5[fCurrentSlice];
1889 fNcs = fN5[fCurrentSlice];
1893 fI = FindClusterIndex(fZmin);
1894 fImax = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
1900 //------------------------------------------------------------------------
1901 Int_t AliITStrackerMI::AliITSlayer::
1902 FindDetectorIndex(Double_t phi, Double_t z) const {
1903 //--------------------------------------------------------------------
1904 //This function finds the detector crossed by the track
1905 //--------------------------------------------------------------------
1907 if (fZOffset<0) // old geometry
1908 dphi = -(phi-fPhiOffset);
1909 else // new geometry
1910 dphi = phi-fPhiOffset;
1913 if (dphi < 0) dphi += 2*TMath::Pi();
1914 else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1915 Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1916 if (np>=fNladders) np-=fNladders;
1917 if (np<0) np+=fNladders;
1920 Double_t dz=fZOffset-z;
1921 Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1922 Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1923 if (nz>=fNdetectors || nz<0) {
1924 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1928 // ad hoc correction for 3rd ladder of SDD inner layer,
1929 // which is reversed (rotated by pi around local y)
1930 // this correction is OK only from AliITSv11Hybrid onwards
1931 if (GetR()>12. && GetR()<20.) { // SDD inner
1932 if(np==2) { // 3rd ladder
1933 Double_t posMod252[3];
1934 AliITSgeomTGeo::GetTranslation(252,posMod252);
1935 // check the Z coordinate of Mod 252: if negative
1936 // (old SDD geometry in AliITSv11Hybrid)
1937 // the swap of numeration whould be applied
1938 if(posMod252[2]<0.){
1939 nz = (fNdetectors-1) - nz;
1943 //printf("ndet %d phi %f z %f np %d nz %d\n",fNdetectors,phi,z,np,nz);
1946 return np*fNdetectors + nz;
1948 //------------------------------------------------------------------------
1949 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1951 //--------------------------------------------------------------------
1952 // This function returns clusters within the "window"
1953 //--------------------------------------------------------------------
1955 if (fCurrentSlice<0) {
1956 Double_t rpi2 = 2.*fR*TMath::Pi();
1957 for (Int_t i=fI; i<fImax; i++) {
1960 if (fYmax<y) y -= rpi2;
1961 if (fYmin>y) y += rpi2;
1962 if (y<fYmin) continue;
1963 if (y>fYmax) continue;
1965 // skip clusters that are in "extended" road but they
1966 // 3sigma error does not touch the original road
1967 if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
1968 if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
1970 if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
1973 return fClusters[i];
1976 for (Int_t i=fI; i<fImax; i++) {
1977 if (fYcs[i]<fYmin) continue;
1978 if (fYcs[i]>fYmax) continue;
1979 if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
1980 ci=fClusterIndexCs[i];
1982 return fClustersCs[i];
1987 //------------------------------------------------------------------------
1988 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1990 //--------------------------------------------------------------------
1991 // This function returns the layer thickness at this point (units X0)
1992 //--------------------------------------------------------------------
1994 x0=AliITSRecoParam::GetX0Air();
1995 if (43<fR&&fR<45) { //SSD2
1998 if (TMath::Abs(y-0.00)>3.40) d+=dd;
1999 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2000 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2001 for (Int_t i=0; i<12; i++) {
2002 if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2003 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2007 if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2008 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2012 if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2013 if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2016 if (37<fR&&fR<41) { //SSD1
2019 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2020 if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2021 if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2022 for (Int_t i=0; i<11; i++) {
2023 if (TMath::Abs(z-3.9*i)<0.15) {
2024 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2028 if (TMath::Abs(z+3.9*i)<0.15) {
2029 if (TMath::Abs(y-0.00)>3.40) d+=dd;
2033 if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2034 if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2037 if (13<fR&&fR<26) { //SDD
2040 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2042 if (TMath::Abs(y-1.80)<0.55) {
2044 for (Int_t j=0; j<20; j++) {
2045 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2046 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2049 if (TMath::Abs(y+1.80)<0.55) {
2051 for (Int_t j=0; j<20; j++) {
2052 if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2053 if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2057 for (Int_t i=0; i<4; i++) {
2058 if (TMath::Abs(z-7.3*i)<0.60) {
2060 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2063 if (TMath::Abs(z+7.3*i)<0.60) {
2065 if (TMath::Abs(y-0.00)>3.30) d+=dd;
2070 if (6<fR&&fR<8) { //SPD2
2071 Double_t dd=0.0063; x0=21.5;
2073 if (TMath::Abs(y-3.08)>0.5) d+=dd;
2074 if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2076 if (3<fR&&fR<5) { //SPD1
2077 Double_t dd=0.0063; x0=21.5;
2079 if (TMath::Abs(y+0.21)>0.6) d+=dd;
2080 if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2085 //------------------------------------------------------------------------
2086 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2088 fRmisal(det.fRmisal),
2090 fSinPhi(det.fSinPhi),
2091 fCosPhi(det.fCosPhi),
2097 fNChips(det.fNChips),
2098 fChipIsBad(det.fChipIsBad)
2102 //------------------------------------------------------------------------
2103 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2104 const AliITSDetTypeRec *detTypeRec)
2106 //--------------------------------------------------------------------
2107 // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2108 //--------------------------------------------------------------------
2110 // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2111 // while in the tracker they start from 0 for each layer
2112 for(Int_t il=0; il<ilayer; il++)
2113 idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2116 if (ilayer==0 || ilayer==1) { // ---------- SPD
2118 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
2120 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
2123 printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2127 // Get calibration from AliITSDetTypeRec
2128 AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2129 calib->SetModuleIndex(idet);
2130 AliITSCalibration *calibSPDdead = 0;
2131 if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2132 if (calib->IsBad() ||
2133 (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2136 // printf("lay %d bad %d\n",ilayer,idet);
2139 // Get segmentation from AliITSDetTypeRec
2140 AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2142 // Read info about bad chips
2143 fNChips = segm->GetMaximumChipIndex()+1;
2144 //printf("ilayer %d detType %d idet %d fNChips %d %d GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2145 if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2146 fChipIsBad = new Bool_t[fNChips];
2147 for (Int_t iCh=0;iCh<fNChips;iCh++) {
2148 fChipIsBad[iCh] = calib->IsChipBad(iCh);
2149 if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2150 //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2155 //------------------------------------------------------------------------
2156 Double_t AliITStrackerMI::GetEffectiveThickness()
2158 //--------------------------------------------------------------------
2159 // Returns the thickness between the current layer and the vertex (units X0)
2160 //--------------------------------------------------------------------
2163 if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2164 if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2165 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2169 Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2170 Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2174 Double_t xn=fgLayers[fI].GetR();
2175 for (Int_t i=0; i<fI; i++) {
2176 Double_t xi=fgLayers[i].GetR();
2177 Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2183 Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2184 d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2187 Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2188 d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2192 //------------------------------------------------------------------------
2193 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2194 //-------------------------------------------------------------------
2195 // This function returns number of clusters within the "window"
2196 //--------------------------------------------------------------------
2198 for (Int_t i=fI; i<fN; i++) {
2199 const AliITSRecPoint *c=fClusters[i];
2200 if (c->GetZ() > fZmax) break;
2201 if (c->IsUsed()) continue;
2202 const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2203 Double_t y=fR*det.GetPhi() + c->GetY();
2205 if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2206 if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2208 if (y<fYmin) continue;
2209 if (y>fYmax) continue;
2214 //------------------------------------------------------------------------
2215 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2216 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2218 //--------------------------------------------------------------------
2219 // This function refits the track "track" at the position "x" using
2220 // the clusters from "clusters"
2221 // If "extra"==kTRUE,
2222 // the clusters from overlapped modules get attached to "track"
2223 // If "planeff"==kTRUE,
2224 // special approach for plane efficiency evaluation is applyed
2225 //--------------------------------------------------------------------
2227 Int_t index[AliITSgeomTGeo::kNLayers];
2229 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2230 Int_t nc=clusters->GetNumberOfClusters();
2231 for (k=0; k<nc; k++) {
2232 Int_t idx=clusters->GetClusterIndex(k);
2233 Int_t ilayer=(idx&0xf0000000)>>28;
2237 return RefitAt(xx,track,index,extra,planeeff); // call the method below
2239 //------------------------------------------------------------------------
2240 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2241 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2243 //--------------------------------------------------------------------
2244 // This function refits the track "track" at the position "x" using
2245 // the clusters from array
2246 // If "extra"==kTRUE,
2247 // the clusters from overlapped modules get attached to "track"
2248 // If "planeff"==kTRUE,
2249 // special approach for plane efficiency evaluation is applyed
2250 //--------------------------------------------------------------------
2251 Int_t index[AliITSgeomTGeo::kNLayers];
2253 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2255 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2256 index[k]=clusters[k];
2259 // special for cosmics and TPC prolonged tracks:
2260 // propagate to the innermost of:
2261 // - innermost layer crossed by the track
2262 // - innermost layer where a cluster was associated to the track
2263 static AliITSRecoParam *repa = NULL;
2265 repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2267 repa = AliITSRecoParam::GetHighFluxParam();
2268 AliWarning("Using default AliITSRecoParam class");
2271 Int_t evsp=repa->GetEventSpecie();
2273 if(track->GetESDtrack()) trStatus=track->GetStatus();
2274 Int_t innermostlayer=0;
2275 if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin)) {
2277 Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2278 for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2279 if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2280 index[innermostlayer] >= 0 ) break;
2283 AliDebug(2,Form(" drphi %f innermost %d",drphi,innermostlayer));
2286 Int_t modstatus=1; // found
2288 Int_t from, to, step;
2289 if (xx > track->GetX()) {
2290 from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2293 from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2296 TString dir = (step>0 ? "outward" : "inward");
2298 for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2299 AliITSlayer &layer=fgLayers[ilayer];
2300 Double_t r=layer.GetR();
2302 if (step<0 && xx>r) break;
2304 // material between SSD and SDD, SDD and SPD
2305 Double_t hI=ilayer-0.5*step;
2306 if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2307 if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2308 if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2309 if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2312 Double_t oldGlobXYZ[3];
2313 if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2315 // continue if we are already beyond this layer
2316 Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2317 if(step>0 && oldGlobR > r) continue; // going outward
2318 if(step<0 && oldGlobR < r) continue; // going inward
2321 if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2323 Int_t idet=layer.FindDetectorIndex(phi,z);
2325 // check if we allow a prolongation without point for large-eta tracks
2326 Int_t skip = CheckSkipLayer(track,ilayer,idet);
2328 modstatus = 4; // out in z
2329 if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2330 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2333 // apply correction for material of the current layer
2334 // add time if going outward
2335 CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2339 if (idet<0) return kFALSE;
2342 const AliITSdetector &det=layer.GetDetector(idet);
2343 // only for ITS-SA tracks refit
2344 if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2346 if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2348 track->SetDetectorIndex(idet);
2349 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2351 Double_t dz,zmin,zmax,dy,ymin,ymax;
2353 const AliITSRecPoint *clAcc=0;
2354 Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2356 Int_t idx=index[ilayer];
2357 if (idx>=0) { // cluster in this layer
2358 modstatus = 6; // no refit
2359 const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2361 if (idet != cl->GetDetectorIndex()) {
2362 idet=cl->GetDetectorIndex();
2363 const AliITSdetector &detc=layer.GetDetector(idet);
2364 if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2365 track->SetDetectorIndex(idet);
2366 if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2368 Int_t cllayer = (idx & 0xf0000000) >> 28;;
2369 Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2373 modstatus = 1; // found
2378 } else { // no cluster in this layer
2380 modstatus = 3; // skipped
2381 // Plane Eff determination:
2382 if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2383 if (IsOKForPlaneEff(track,clusters,ilayer)) // only adequate track for plane eff. evaluation
2384 UseTrackForPlaneEff(track,ilayer);
2387 modstatus = 5; // no cls in road
2389 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2390 dz = 0.5*(zmax-zmin);
2391 dy = 0.5*(ymax-ymin);
2392 Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2393 if (dead==1) modstatus = 7; // holes in z in SPD
2394 if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2399 if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2400 track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2402 track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2405 if (extra && clAcc) { // search for extra clusters in overlapped modules
2406 AliITStrackV2 tmp(*track);
2407 if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2408 layer.SelectClusters(zmin,zmax,ymin,ymax);
2410 const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2412 maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2413 Double_t tolerance=0.1;
2414 while ((clExtra=layer.GetNextCluster(ci))!=0) {
2415 // only clusters in another module! (overlaps)
2416 idetExtra = clExtra->GetDetectorIndex();
2417 if (idet == idetExtra) continue;
2419 const AliITSdetector &detx=layer.GetDetector(idetExtra);
2421 if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2422 if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2423 if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2424 if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2426 Double_t chi2=tmp.GetPredictedChi2(clExtra);
2427 if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2430 track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2431 track->SetExtraModule(ilayer,idetExtra);
2433 } // end search for extra clusters in overlapped modules
2435 // Correct for material of the current layer
2437 // add time if going outward
2438 if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2439 track->SetCheckInvariant(kTRUE);
2440 } // end loop on layers
2442 if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2446 //------------------------------------------------------------------------
2447 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2450 // calculate normalized chi2
2451 // return NormalizedChi2(track,0);
2454 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2455 // track->fdEdxMismatch=0;
2456 Float_t dedxmismatch =0;
2457 Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2459 for (Int_t i = 0;i<6;i++){
2460 if (track->GetClIndex(i)>=0){
2461 Float_t cerry, cerrz;
2462 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2464 { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2467 Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2468 if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2469 Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2471 cchi2+=(0.5-ratio)*10.;
2472 //track->fdEdxMismatch+=(0.5-ratio)*10.;
2473 dedxmismatch+=(0.5-ratio)*10.;
2477 AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2478 Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2479 if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2480 if (i<2) chi2+=2*cl->GetDeltaProbability();
2486 if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2487 track->SetdEdxMismatch(dedxmismatch);
2491 for (Int_t i = 0;i<4;i++){
2492 if (track->GetClIndex(i)>=0){
2493 Float_t cerry, cerrz;
2494 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2495 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2498 chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2499 chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2503 for (Int_t i = 4;i<6;i++){
2504 if (track->GetClIndex(i)>=0){
2505 Float_t cerry, cerrz;
2506 if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2507 else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2510 Float_t cerryb, cerrzb;
2511 if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2512 else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2515 chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2516 chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2521 if (track->GetESDtrack()->GetTPCsignal()>85){
2522 Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2524 chi2+=(0.5-ratio)*5.;
2527 chi2+=(ratio-2.0)*3;
2531 Double_t match = TMath::Sqrt(track->GetChi22());
2532 if (track->GetConstrain()) match/=track->GetNumberOfClusters();
2533 if (!track->GetConstrain()) {
2534 if (track->GetNumberOfClusters()>2) {
2535 match/=track->GetNumberOfClusters()-2.;
2540 if (match<0) match=0;
2542 // penalty factor for missing points (NDeadZone>0), but no penalty
2543 // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2544 // or there is a dead from OCDB)
2545 Float_t deadzonefactor = 0.;
2546 if (track->GetNDeadZone()>0.) {
2547 Int_t sumDeadZoneProbability=0;
2548 for(Int_t ilay=0;ilay<6;ilay++) {
2549 if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2551 Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2552 if(nDeadZoneWithProbNot1>0) {
2553 Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2554 AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2555 deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2557 deadZoneProbability = TMath::Min(deadZoneProbability,one);
2558 deadzonefactor = 3.*(1.1-deadZoneProbability);
2562 Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2563 (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2564 1./(1.+track->GetNSkipped()));
2565 AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2566 AliDebug(2,Form("NormChi2 %f cls %d\n",normchi2,track->GetNumberOfClusters()));
2569 //------------------------------------------------------------------------
2570 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2573 // return matching chi2 between two tracks
2574 Double_t largeChi2=1000.;
2576 AliITStrackMI track3(*track2);
2577 if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2579 vec(0,0)=track1->GetY() - track3.GetY();
2580 vec(1,0)=track1->GetZ() - track3.GetZ();
2581 vec(2,0)=track1->GetSnp() - track3.GetSnp();
2582 vec(3,0)=track1->GetTgl() - track3.GetTgl();
2583 vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2586 cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2587 cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2588 cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2589 cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2590 cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2592 cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2593 cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2594 cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2595 cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2597 cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2598 cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2599 cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2601 cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2602 cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2604 cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2607 TMatrixD vec2(cov,TMatrixD::kMult,vec);
2608 TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2611 //------------------------------------------------------------------------
2612 Double_t AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2615 // return probability that given point (characterized by z position and error)
2616 // is in SPD dead zone
2617 // This method assumes that fSPDdetzcentre is ordered from -z to +z
2619 Double_t probability = 0.;
2620 Double_t nearestz = 0.,distz=0.;
2621 Int_t nearestzone = -1;
2622 Double_t mindistz = 1000.;
2624 // find closest dead zone
2625 for (Int_t i=0; i<3; i++) {
2626 distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2627 if (distz<mindistz) {
2629 nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2634 // too far from dead zone
2635 if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2638 Double_t zmin, zmax;
2639 if (nearestzone==0) { // dead zone at z = -7
2640 zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2641 zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2642 } else if (nearestzone==1) { // dead zone at z = 0
2643 zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2644 zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2645 } else if (nearestzone==2) { // dead zone at z = +7
2646 zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2647 zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2652 // probability that the true z is in the range [zmin,zmax] (i.e. inside
2654 probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2655 AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2656 AliDebug(2,Form("zpos %f +- %f nearestzone %d zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2659 //------------------------------------------------------------------------
2660 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2663 // calculate normalized chi2
2665 Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2667 for (Int_t i = 0;i<6;i++){
2668 if (TMath::Abs(track->GetDy(i))>0){
2669 chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2670 chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2673 else{chi2[i]=10000;}
2676 TMath::Sort(6,chi2,index,kFALSE);
2677 Float_t max = float(ncl)*fac-1.;
2678 Float_t sumchi=0, sumweight=0;
2679 for (Int_t i=0;i<max+1;i++){
2680 Float_t weight = (i<max)?1.:(max+1.-i);
2681 sumchi+=weight*chi2[index[i]];
2684 Double_t normchi2 = sumchi/sumweight;
2687 //------------------------------------------------------------------------
2688 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2691 // calculate normalized chi2
2692 // if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2695 for (Int_t i=0;i<6;i++){
2696 if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2697 Double_t sy1 = forwardtrack->GetSigmaY(i);
2698 Double_t sz1 = forwardtrack->GetSigmaZ(i);
2699 Double_t sy2 = backtrack->GetSigmaY(i);
2700 Double_t sz2 = backtrack->GetSigmaZ(i);
2701 if (i<2){ sy2=1000.;sz2=1000;}
2703 Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2704 Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2706 Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2707 Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2709 res+= nz0*nz0+ny0*ny0;
2712 if (npoints>1) return
2713 TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2714 //2*forwardtrack->fNUsed+
2715 res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2716 1./(1.+forwardtrack->GetNSkipped()));
2719 //------------------------------------------------------------------------
2720 Float_t *AliITStrackerMI::GetWeight(Int_t index) {
2721 //--------------------------------------------------------------------
2722 // Return pointer to a given cluster
2723 //--------------------------------------------------------------------
2724 Int_t l=(index & 0xf0000000) >> 28;
2725 Int_t c=(index & 0x0fffffff) >> 00;
2726 return fgLayers[l].GetWeight(c);
2728 //------------------------------------------------------------------------
2729 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2731 //---------------------------------------------
2732 // register track to the list
2734 if (track->GetESDtrack()->GetKinkIndex(0)!=0) return; //don't register kink tracks
2737 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2738 Int_t index = track->GetClusterIndex(icluster);
2739 Int_t l=(index & 0xf0000000) >> 28;
2740 Int_t c=(index & 0x0fffffff) >> 00;
2741 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2742 for (Int_t itrack=0;itrack<4;itrack++){
2743 if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2744 fgLayers[l].SetClusterTracks(itrack,c,id);
2750 //------------------------------------------------------------------------
2751 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2753 //---------------------------------------------
2754 // unregister track from the list
2755 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2756 Int_t index = track->GetClusterIndex(icluster);
2757 Int_t l=(index & 0xf0000000) >> 28;
2758 Int_t c=(index & 0x0fffffff) >> 00;
2759 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2760 for (Int_t itrack=0;itrack<4;itrack++){
2761 if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2762 fgLayers[l].SetClusterTracks(itrack,c,-1);
2767 //------------------------------------------------------------------------
2768 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2770 //-------------------------------------------------------------
2771 //get number of shared clusters
2772 //-------------------------------------------------------------
2774 for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2775 // mean number of clusters
2776 Float_t *ny = GetNy(id), *nz = GetNz(id);
2779 for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2780 Int_t index = track->GetClusterIndex(icluster);
2781 Int_t l=(index & 0xf0000000) >> 28;
2782 Int_t c=(index & 0x0fffffff) >> 00;
2783 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2784 // if (ny[l]<1.e-13){
2785 // printf("problem\n");
2787 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2791 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2792 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2793 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2795 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2798 deltan = (cl->GetNz()-nz[l]);
2800 if (deltan>2.0) continue; // extended - highly probable shared cluster
2801 weight = 2./TMath::Max(3.+deltan,2.);
2803 for (Int_t itrack=0;itrack<4;itrack++){
2804 if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
2806 clist[l] = (AliITSRecPoint*)GetCluster(index);
2812 track->SetNUsed(shared);
2815 //------------------------------------------------------------------------
2816 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
2819 // find first shared track
2821 // mean number of clusters
2822 Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
2824 for (Int_t i=0;i<6;i++) overlist[i]=-1;
2825 Int_t sharedtrack=100000;
2826 Int_t tracks[24],trackindex=0;
2827 for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
2829 for (Int_t icluster=0;icluster<6;icluster++){
2830 if (clusterlist[icluster]<0) continue;
2831 Int_t index = clusterlist[icluster];
2832 Int_t l=(index & 0xf0000000) >> 28;
2833 Int_t c=(index & 0x0fffffff) >> 00;
2834 // if (ny[l]<1.e-13){
2835 // printf("problem\n");
2837 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2838 //if (l>3) continue;
2839 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2842 if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
2843 if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
2844 if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
2846 deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
2849 deltan = (cl->GetNz()-nz[l]);
2851 if (deltan>2.0) continue; // extended - highly probable shared cluster
2853 for (Int_t itrack=3;itrack>=0;itrack--){
2854 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2855 if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
2856 tracks[trackindex] = fgLayers[l].GetClusterTracks(itrack,c);
2861 if (trackindex==0) return -1;
2863 sharedtrack = tracks[0];
2865 if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
2868 Int_t tracks2[24], cluster[24];
2869 for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
2872 for (Int_t i=0;i<trackindex;i++){
2873 if (tracks[i]<0) continue;
2874 tracks2[index] = tracks[i];
2876 for (Int_t j=i+1;j<trackindex;j++){
2877 if (tracks[j]<0) continue;
2878 if (tracks[j]==tracks[i]){
2886 for (Int_t i=0;i<index;i++){
2887 if (cluster[index]>max) {
2888 sharedtrack=tracks2[index];
2895 if (sharedtrack>=100000) return -1;
2897 // make list of overlaps
2899 for (Int_t icluster=0;icluster<6;icluster++){
2900 if (clusterlist[icluster]<0) continue;
2901 Int_t index = clusterlist[icluster];
2902 Int_t l=(index & 0xf0000000) >> 28;
2903 Int_t c=(index & 0x0fffffff) >> 00;
2904 if (c>fgLayers[l].GetNumberOfClusters()) continue;
2905 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
2907 if (cl->GetNy()>2) continue;
2908 if (cl->GetNz()>2) continue;
2911 if (cl->GetNy()>3) continue;
2912 if (cl->GetNz()>3) continue;
2915 for (Int_t itrack=3;itrack>=0;itrack--){
2916 if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
2917 if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
2925 //------------------------------------------------------------------------
2926 AliITStrackMI * AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
2928 // try to find track hypothesys without conflicts
2929 // with minimal chi2;
2930 TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
2931 Int_t entries1 = arr1->GetEntriesFast();
2932 TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
2933 if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
2934 Int_t entries2 = arr2->GetEntriesFast();
2935 if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
2937 AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
2938 AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
2939 if (track10->Pt()>0.5+track20->Pt()) return track10;
2941 for (Int_t itrack=0;itrack<entries1;itrack++){
2942 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2943 UnRegisterClusterTracks(track,trackID1);
2946 for (Int_t itrack=0;itrack<entries2;itrack++){
2947 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2948 UnRegisterClusterTracks(track,trackID2);
2952 Float_t maxconflicts=6;
2953 Double_t maxchi2 =1000.;
2955 // get weight of hypothesys - using best hypothesy
2958 Int_t list1[6],list2[6];
2959 AliITSRecPoint *clist1[6], *clist2[6] ;
2960 RegisterClusterTracks(track10,trackID1);
2961 RegisterClusterTracks(track20,trackID2);
2962 Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
2963 Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
2964 UnRegisterClusterTracks(track10,trackID1);
2965 UnRegisterClusterTracks(track20,trackID2);
2968 Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
2969 Float_t nerry[6],nerrz[6];
2970 Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
2971 Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
2972 for (Int_t i=0;i<6;i++){
2973 if ( (erry1[i]>0) && (erry2[i]>0)) {
2974 nerry[i] = TMath::Min(erry1[i],erry2[i]);
2975 nerrz[i] = TMath::Min(errz1[i],errz2[i]);
2977 nerry[i] = TMath::Max(erry1[i],erry2[i]);
2978 nerrz[i] = TMath::Max(errz1[i],errz2[i]);
2980 if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
2981 chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
2982 chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
2985 if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
2986 chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
2987 chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
2995 Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
2996 Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
2997 Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
2998 Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3000 w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3001 +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3002 +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3004 w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3005 s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3006 +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3009 Double_t sumw = w1+w2;
3013 w1 = (d2+0.5)/(d1+d2+1.);
3014 w2 = (d1+0.5)/(d1+d2+1.);
3016 // Float_t maxmax = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3017 //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3019 // get pair of "best" hypothesys
3021 Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3022 Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3024 for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3025 AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3026 //if (track1->fFakeRatio>0) continue;
3027 RegisterClusterTracks(track1,trackID1);
3028 for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3029 AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3031 // Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3032 //if (track2->fFakeRatio>0) continue;
3034 RegisterClusterTracks(track2,trackID2);
3035 Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3036 Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3037 UnRegisterClusterTracks(track2,trackID2);
3039 if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3040 if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3041 if (nskipped>0.5) continue;
3043 //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3044 if (conflict1+1<cconflict1) continue;
3045 if (conflict2+1<cconflict2) continue;
3049 for (Int_t i=0;i<6;i++){
3051 Float_t c1 =0.; // conflict coeficients
3053 if (clist1[i]&&clist2[i]){
3056 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3059 deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3061 c1 = 2./TMath::Max(3.+deltan,2.);
3062 c2 = 2./TMath::Max(3.+deltan,2.);
3068 deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3071 deltan = (clist1[i]->GetNz()-nz1[i]);
3073 c1 = 2./TMath::Max(3.+deltan,2.);
3080 deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3083 deltan = (clist2[i]->GetNz()-nz2[i]);
3085 c2 = 2./TMath::Max(3.+deltan,2.);
3091 if (TMath::Abs(track1->GetDy(i))>0.) {
3092 chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3093 (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3094 //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3095 // (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3097 if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3100 if (TMath::Abs(track2->GetDy(i))>0.) {
3101 chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3102 (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3103 //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3104 // (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3107 if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3109 sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3110 if (chi21>0) sum+=w1;
3111 if (chi22>0) sum+=w2;
3114 Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3115 if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3116 Double_t normchi2 = 2*conflict+sumchi2/sum;
3117 if ( normchi2 <maxchi2 ){
3120 maxconflicts = conflict;
3124 UnRegisterClusterTracks(track1,trackID1);
3127 // if (maxconflicts<4 && maxchi2<th0){
3128 if (maxchi2<th0*2.){
3129 Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3130 AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3131 track1->SetChi2MIP(5,maxconflicts);
3132 track1->SetChi2MIP(6,maxchi2);
3133 track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3134 // track1->UpdateESDtrack(AliESDtrack::kITSin);
3135 track1->SetChi2MIP(8,index1);
3136 fBestTrackIndex[trackID1] =index1;
3137 UpdateESDtrack(track1, AliESDtrack::kITSin);
3139 else if (track10->GetChi2MIP(0)<th1){
3140 track10->SetChi2MIP(5,maxconflicts);
3141 track10->SetChi2MIP(6,maxchi2);
3142 // track10->UpdateESDtrack(AliESDtrack::kITSin);
3143 UpdateESDtrack(track10,AliESDtrack::kITSin);
3146 for (Int_t itrack=0;itrack<entries1;itrack++){
3147 AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3148 UnRegisterClusterTracks(track,trackID1);
3151 for (Int_t itrack=0;itrack<entries2;itrack++){
3152 AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3153 UnRegisterClusterTracks(track,trackID2);
3156 if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3157 &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3158 // if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3159 // &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3160 RegisterClusterTracks(track10,trackID1);
3162 if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3163 &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3164 //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3165 // &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3166 RegisterClusterTracks(track20,trackID2);
3171 //------------------------------------------------------------------------
3172 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3173 //--------------------------------------------------------------------
3174 // This function marks clusters assigned to the track
3175 //--------------------------------------------------------------------
3176 AliTracker::UseClusters(t,from);
3178 AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3179 //if (c->GetQ()>2) c->Use();
3180 if (c->GetSigmaZ2()>0.1) c->Use();
3181 c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3182 //if (c->GetQ()>2) c->Use();
3183 if (c->GetSigmaZ2()>0.1) c->Use();
3186 //------------------------------------------------------------------------
3187 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3189 //------------------------------------------------------------------
3190 // add track to the list of hypothesys
3191 //------------------------------------------------------------------
3193 if (esdindex>=fTrackHypothesys.GetEntriesFast())
3194 fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
3196 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3198 array = new TObjArray(10);
3199 fTrackHypothesys.AddAt(array,esdindex);
3201 array->AddLast(track);
3203 //------------------------------------------------------------------------
3204 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3206 //-------------------------------------------------------------------
3207 // compress array of track hypothesys
3208 // keep only maxsize best hypothesys
3209 //-------------------------------------------------------------------
3210 if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3211 if (! (fTrackHypothesys.At(esdindex)) ) return;
3212 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3213 Int_t entries = array->GetEntriesFast();
3215 //- find preliminary besttrack as a reference
3216 Float_t minchi2=10000;
3218 AliITStrackMI * besttrack=0;
3219 for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3220 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3221 if (!track) continue;
3222 Float_t chi2 = NormalizedChi2(track,0);
3224 Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3225 track->SetLabel(tpcLabel);
3227 track->SetFakeRatio(1.);
3228 CookLabel(track,0.); //For comparison only
3230 //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3231 if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3232 if (track->GetNumberOfClusters()<maxn) continue;
3233 maxn = track->GetNumberOfClusters();
3240 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3241 delete array->RemoveAt(itrack);
3245 if (!besttrack) return;
3248 //take errors of best track as a reference
3249 Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3250 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3251 for (Int_t j=0;j<6;j++) {
3252 if (besttrack->GetClIndex(j)>=0){
3253 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3254 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3255 ny[j] = besttrack->GetNy(j);
3256 nz[j] = besttrack->GetNz(j);
3260 // calculate normalized chi2
3262 Float_t * chi2 = new Float_t[entries];
3263 Int_t * index = new Int_t[entries];
3264 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3265 for (Int_t itrack=0;itrack<entries;itrack++){
3266 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3268 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3269 track->SetChi2MIP(0,GetNormalizedChi2(track, mode));
3270 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3271 chi2[itrack] = track->GetChi2MIP(0);
3273 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3274 delete array->RemoveAt(itrack);
3280 TMath::Sort(entries,chi2,index,kFALSE);
3281 besttrack = (AliITStrackMI*)array->At(index[0]);
3282 if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3283 if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3284 for (Int_t j=0;j<6;j++){
3285 if (besttrack->GetClIndex(j)>=0){
3286 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3287 errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3288 ny[j] = besttrack->GetNy(j);
3289 nz[j] = besttrack->GetNz(j);
3294 // calculate one more time with updated normalized errors
3295 for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3296 for (Int_t itrack=0;itrack<entries;itrack++){
3297 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3299 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3300 AliDebug(2,Form("track %d ncls %d\n",itrack,track->GetNumberOfClusters()));
3301 if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0))
3302 chi2[itrack] = track->GetChi2MIP(0)-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3305 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3306 delete array->RemoveAt(itrack);
3311 entries = array->GetEntriesFast();
3315 TObjArray * newarray = new TObjArray();
3316 TMath::Sort(entries,chi2,index,kFALSE);
3317 besttrack = (AliITStrackMI*)array->At(index[0]);
3319 AliDebug(2,Form("ncls best track %d %f %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3321 for (Int_t j=0;j<6;j++){
3322 if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3323 erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3324 errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3325 ny[j] = besttrack->GetNy(j);
3326 nz[j] = besttrack->GetNz(j);
3329 besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3330 minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3331 Float_t minn = besttrack->GetNumberOfClusters()-3;
3333 for (Int_t i=0;i<entries;i++){
3334 AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3335 if (!track) continue;
3336 if (accepted>maxcut) break;
3337 track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3338 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3339 if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3340 delete array->RemoveAt(index[i]);
3344 Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3345 if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3346 if (!shortbest) accepted++;
3348 newarray->AddLast(array->RemoveAt(index[i]));
3349 for (Int_t j=0;j<6;j++){
3351 erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3352 errz[j] = track->GetSigmaZ(j); errz[j] = track->GetSigmaZ(j+6);
3353 ny[j] = track->GetNy(j);
3354 nz[j] = track->GetNz(j);
3359 delete array->RemoveAt(index[i]);
3363 delete fTrackHypothesys.RemoveAt(esdindex);
3364 fTrackHypothesys.AddAt(newarray,esdindex);
3368 delete fTrackHypothesys.RemoveAt(esdindex);
3374 //------------------------------------------------------------------------
3375 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3377 //-------------------------------------------------------------
3378 // try to find best hypothesy
3379 // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3380 //-------------------------------------------------------------
3381 if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3382 TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3383 if (!array) return 0;
3384 Int_t entries = array->GetEntriesFast();
3385 if (!entries) return 0;
3386 Float_t minchi2 = 100000;
3387 AliITStrackMI * besttrack=0;
3389 AliITStrackMI * backtrack = new AliITStrackMI(*original);
3390 AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3391 Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3392 Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3394 for (Int_t i=0;i<entries;i++){
3395 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3396 if (!track) continue;
3397 Float_t sigmarfi,sigmaz;
3398 GetDCASigma(track,sigmarfi,sigmaz);
3399 track->SetDnorm(0,sigmarfi);
3400 track->SetDnorm(1,sigmaz);
3402 track->SetChi2MIP(1,1000000);
3403 track->SetChi2MIP(2,1000000);
3404 track->SetChi2MIP(3,1000000);
3407 backtrack = new(backtrack) AliITStrackMI(*track);
3408 if (track->GetConstrain()) {
3409 if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3410 if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3411 if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;
3413 backtrack->ResetCovariance(10.);
3415 backtrack->ResetCovariance(10.);
3417 backtrack->ResetClusters();
3419 Double_t x = original->GetX();
3420 if (!RefitAt(x,backtrack,track)) continue;
3422 track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3423 //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3424 if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.) continue;
3425 track->SetChi22(GetMatchingChi2(backtrack,original));
3427 if ((track->GetConstrain()) && track->GetChi22()>90.) continue;
3428 if ((!track->GetConstrain()) && track->GetChi22()>30.) continue;
3429 if ( track->GetChi22()/track->GetNumberOfClusters()>11.) continue;
3432 if (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)) continue;
3434 //forward track - without constraint
3435 forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3436 forwardtrack->ResetClusters();
3438 RefitAt(x,forwardtrack,track);
3439 track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3440 if (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0) continue;
3441 if (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)) continue;
3443 //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3444 //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3445 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP()); //I.B.
3446 forwardtrack->SetD(0,track->GetD(0));
3447 forwardtrack->SetD(1,track->GetD(1));
3450 AliITSRecPoint* clist[6];
3451 track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3452 if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3455 track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3456 if ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3457 if ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3458 track->SetChi2MIP(3,1000);
3461 Double_t chi2 = track->GetChi2MIP(0)+track->GetNUsed();
3463 for (Int_t ichi=0;ichi<5;ichi++){
3464 forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3466 if (chi2 < minchi2){
3467 //besttrack = new AliITStrackMI(*forwardtrack);
3469 besttrack->SetLabel(track->GetLabel());
3470 besttrack->SetFakeRatio(track->GetFakeRatio());
3472 //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3473 //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3474 forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP()); //I.B.
3478 delete forwardtrack;
3480 for (Int_t i=0;i<entries;i++){
3481 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3483 if (!track) continue;
3485 if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3486 (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
3487 track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.){
3488 if (track->GetConstrain() || track->GetNumberOfClusters()>5){ //keep best short tracks - without vertex constrain
3489 delete array->RemoveAt(i);
3499 SortTrackHypothesys(esdindex,checkmax,1);
3501 array = (TObjArray*) fTrackHypothesys.At(esdindex);
3502 if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3503 besttrack = (AliITStrackMI*)array->At(0);
3504 if (!besttrack) return 0;
3505 besttrack->SetChi2MIP(8,0);
3506 fBestTrackIndex[esdindex]=0;
3507 entries = array->GetEntriesFast();
3508 AliITStrackMI *longtrack =0;
3510 Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3511 for (Int_t itrack=entries-1;itrack>0;itrack--) {
3512 AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3513 if (!track->GetConstrain()) continue;
3514 if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3515 if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3516 if (track->GetChi2MIP(0)>4.) continue;
3517 minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3520 //if (longtrack) besttrack=longtrack;
3523 AliITSRecPoint * clist[6];
3524 Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3525 if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3526 &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3527 RegisterClusterTracks(besttrack,esdindex);
3534 Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3535 if (sharedtrack>=0){
3537 besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
3539 shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3545 if (shared>2.5) return 0;
3546 if (shared>1.0) return besttrack;
3548 // Don't sign clusters if not gold track
3550 if (!besttrack->IsGoldPrimary()) return besttrack;
3551 if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack; //track belong to kink
3553 if (fConstraint[fPass]){
3557 Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3558 for (Int_t i=0;i<6;i++){
3559 Int_t index = besttrack->GetClIndex(i);
3560 if (index<0) continue;
3561 Int_t ilayer = (index & 0xf0000000) >> 28;
3562 if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3563 AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3565 if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3566 if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3567 if ( c->GetNz()> nz[i]+0.7) continue; //shared track
3568 if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3569 if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3570 //if ( c->GetNy()> ny[i]+0.7) continue; //shared track
3572 Bool_t cansign = kTRUE;
3573 for (Int_t itrack=0;itrack<entries; itrack++){
3574 AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3575 if (!track) continue;
3576 if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3577 if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3583 if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3584 if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3585 if (!c->IsUsed()) c->Use();
3591 //------------------------------------------------------------------------
3592 void AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3595 // get "best" hypothesys
3598 Int_t nentries = itsTracks.GetEntriesFast();
3599 for (Int_t i=0;i<nentries;i++){
3600 AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3601 if (!track) continue;
3602 TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3603 if (!array) continue;
3604 if (array->GetEntriesFast()<=0) continue;
3606 AliITStrackMI* longtrack=0;
3608 Float_t maxchi2=1000;
3609 for (Int_t j=0;j<array->GetEntriesFast();j++){
3610 AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3611 if (!trackHyp) continue;
3612 if (trackHyp->GetGoldV0()) {
3613 longtrack = trackHyp; //gold V0 track taken
3616 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3617 Float_t chi2 = trackHyp->GetChi2MIP(0);
3619 if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3621 if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = trackHyp->GetChi2MIP(0);
3623 if (chi2 > maxchi2) continue;
3624 minn= trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3631 AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3632 if (!longtrack) {longtrack = besttrack;}
3633 else besttrack= longtrack;
3637 AliITSRecPoint * clist[6];
3638 Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3640 track->SetNUsed(shared);
3641 track->SetNSkipped(besttrack->GetNSkipped());
3642 track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3644 if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3648 Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3649 //if (sharedtrack==-1) sharedtrack=0;
3650 if (sharedtrack>=0) {
3651 besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
3654 if (besttrack&&fAfterV0) {
3655 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3657 if (besttrack&&fConstraint[fPass])
3658 UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3659 if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3660 if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3661 TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3667 //------------------------------------------------------------------------
3668 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3669 //--------------------------------------------------------------------
3670 //This function "cooks" a track label. If label<0, this track is fake.
3671 //--------------------------------------------------------------------
3674 if (track->GetESDtrack()){
3675 tpcLabel = track->GetESDtrack()->GetTPCLabel();
3676 ULong_t trStatus=track->GetESDtrack()->GetStatus();
3677 if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
3679 track->SetChi2MIP(9,0);
3681 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3682 Int_t cindex = track->GetClusterIndex(i);
3683 Int_t l=(cindex & 0xf0000000) >> 28;
3684 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3686 for (Int_t ind=0;ind<3;ind++){
3687 if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
3688 //AliDebug(2,Form("icl %d ilab %d lab %d",i,ind,cl->GetLabel(ind)));
3690 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3693 Int_t nclusters = track->GetNumberOfClusters();
3694 if (nclusters > 0) //PH Some tracks don't have any cluster
3695 track->SetFakeRatio(double(nwrong)/double(nclusters));
3696 if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
3697 track->SetLabel(-tpcLabel);
3699 track->SetLabel(tpcLabel);
3701 AliDebug(2,Form(" nls %d wrong %d label %d tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
3704 //------------------------------------------------------------------------
3705 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
3707 // Fill the dE/dx in this track
3709 track->SetChi2MIP(9,0);
3710 for (Int_t i=0;i<track->GetNumberOfClusters();i++){
3711 Int_t cindex = track->GetClusterIndex(i);
3712 Int_t l=(cindex & 0xf0000000) >> 28;
3713 AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
3714 Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
3716 for (Int_t ind=0;ind<3;ind++){
3717 if (cl->GetLabel(ind)==lab) isWrong=0;
3719 track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
3723 track->CookdEdx(low,up);
3725 //------------------------------------------------------------------------
3726 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
3728 // Create some arrays
3730 if (fCoefficients) delete []fCoefficients;
3731 fCoefficients = new Float_t[ntracks*48];
3732 for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
3734 //------------------------------------------------------------------------
3735 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
3738 // Compute predicted chi2
3740 // Take into account the mis-alignment (bring track to cluster plane)
3741 Double_t xTrOrig=track->GetX();
3742 if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
3743 Float_t erry,errz,covyz;
3744 Float_t theta = track->GetTgl();
3745 Float_t phi = track->GetSnp();
3746 phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
3747 AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
3748 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3749 AliDebug(2,Form(" chi2: tr-cl %f %f tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
3750 Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
3751 // Bring the track back to detector plane in ideal geometry
3752 // [mis-alignment will be accounted for in UpdateMI()]
3753 if (!track->Propagate(xTrOrig)) return 1000.;
3755 AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
3756 Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
3758 chi2+=0.5*TMath::Min(delta/2,2.);
3759 chi2+=2.*cluster->GetDeltaProbability();
3762 track->SetNy(layer,ny);
3763 track->SetNz(layer,nz);
3764 track->SetSigmaY(layer,erry);
3765 track->SetSigmaZ(layer, errz);
3766 track->SetSigmaYZ(layer,covyz);
3767 //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
3768 track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
3772 //------------------------------------------------------------------------
3773 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
3778 Int_t layer = (index & 0xf0000000) >> 28;
3779 track->SetClIndex(layer, index);
3780 if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
3781 if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
3782 chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
3783 track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
3787 if (TMath::Abs(cl->GetQ())<1.e-13) return 0; // ingore the "virtual" clusters
3790 // Take into account the mis-alignment (bring track to cluster plane)
3791 Double_t xTrOrig=track->GetX();
3792 Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
3793 AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
3794 AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
3795 AliDebug(2,Form(" xtr %f xcl %f",track->GetX(),cl->GetX()));
3797 if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
3800 c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
3801 c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
3802 c.SetSigmaYZ(track->GetSigmaYZ(layer));
3805 Int_t updated = track->UpdateMI(&c,chi2,index);
3806 // Bring the track back to detector plane in ideal geometry
3807 if (!track->Propagate(xTrOrig)) return 0;
3809 if(!updated) AliDebug(2,"update failed");
3813 //------------------------------------------------------------------------
3814 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
3817 //DCA sigmas parameterization
3818 //to be paramterized using external parameters in future
3821 sigmarfi = 0.0040+1.4 *TMath::Abs(track->GetC())+332.*track->GetC()*track->GetC();
3822 sigmaz = 0.0110+4.37*TMath::Abs(track->GetC());
3824 //------------------------------------------------------------------------
3825 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
3828 // Clusters from delta electrons?
3830 Int_t entries = clusterArray->GetEntriesFast();
3831 if (entries<4) return;
3832 AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
3833 Int_t layer = cluster->GetLayer();
3834 if (layer>1) return;
3836 Int_t ncandidates=0;
3837 Float_t r = (layer>0)? 7:4;
3839 for (Int_t i=0;i<entries;i++){
3840 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
3841 Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
3842 if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
3843 index[ncandidates] = i; //candidate to belong to delta electron track
3845 if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3846 cl0->SetDeltaProbability(1);
3852 for (Int_t i=0;i<ncandidates;i++){
3853 AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
3854 if (cl0->GetDeltaProbability()>0.8) continue;
3857 Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3858 sumy=sumz=sumy2=sumyz=sumw=0.0;
3859 for (Int_t j=0;j<ncandidates;j++){
3861 AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
3863 Float_t dz = cl0->GetZ()-cl1->GetZ();
3864 Float_t dy = cl0->GetY()-cl1->GetY();
3865 if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3866 Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3867 y[ncl] = cl1->GetY();
3868 z[ncl] = cl1->GetZ();
3869 sumy+= y[ncl]*weight;
3870 sumz+= z[ncl]*weight;
3871 sumy2+=y[ncl]*y[ncl]*weight;
3872 sumyz+=y[ncl]*z[ncl]*weight;
3877 if (ncl<4) continue;
3878 Float_t det = sumw*sumy2 - sumy*sumy;
3880 if (TMath::Abs(det)>0.01){
3881 Float_t z0 = (sumy2*sumz - sumy*sumyz)/det;
3882 Float_t k = (sumyz*sumw - sumy*sumz)/det;
3883 delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3886 Float_t z0 = sumyz/sumy;
3887 delta = TMath::Abs(cl0->GetZ()-z0);
3890 cl0->SetDeltaProbability(1-20.*delta);
3894 //------------------------------------------------------------------------
3895 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3900 track->UpdateESDtrack(flags);
3901 AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
3902 if (oldtrack) delete oldtrack;
3903 track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
3904 // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
3905 // printf("Problem\n");
3908 //------------------------------------------------------------------------
3909 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
3911 // Get nearest upper layer close to the point xr.
3912 // rough approximation
3914 const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
3915 Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
3917 for (Int_t i=0;i<6;i++){
3918 if (radius<kRadiuses[i]){
3925 //------------------------------------------------------------------------
3926 void AliITStrackerMI::BuildMaterialLUT(TString material) {
3927 //--------------------------------------------------------------------
3928 // Fill a look-up table with mean material
3929 //--------------------------------------------------------------------
3933 Double_t point1[3],point2[3];
3934 Double_t phi,cosphi,sinphi,z;
3935 // 0-5 layers, 6 pipe, 7-8 shields
3936 Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
3937 Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
3939 Int_t ifirst=0,ilast=0;
3940 if(material.Contains("Pipe")) {
3942 } else if(material.Contains("Shields")) {
3944 } else if(material.Contains("Layers")) {
3947 Error("BuildMaterialLUT","Wrong layer name\n");
3950 for(Int_t imat=ifirst; imat<=ilast; imat++) {
3951 Double_t param[5]={0.,0.,0.,0.,0.};
3952 for (Int_t i=0; i<n; i++) {
3953 phi = 2.*TMath::Pi()*gRandom->Rndm();
3954 cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
3955 z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
3956 point1[0] = rmin[imat]*cosphi;
3957 point1[1] = rmin[imat]*sinphi;
3959 point2[0] = rmax[imat]*cosphi;
3960 point2[1] = rmax[imat]*sinphi;
3962 AliTracker::MeanMaterialBudget(point1,point2,mparam);
3963 for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
3965 for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
3967 fxOverX0Layer[imat] = param[1];
3968 fxTimesRhoLayer[imat] = param[0]*param[4];
3969 } else if(imat==6) {
3970 fxOverX0Pipe = param[1];
3971 fxTimesRhoPipe = param[0]*param[4];
3972 } else if(imat==7) {
3973 fxOverX0Shield[0] = param[1];
3974 fxTimesRhoShield[0] = param[0]*param[4];
3975 } else if(imat==8) {
3976 fxOverX0Shield[1] = param[1];
3977 fxTimesRhoShield[1] = param[0]*param[4];
3981 printf("%s\n",material.Data());
3982 printf("%f %f\n",fxOverX0Pipe,fxTimesRhoPipe);
3983 printf("%f %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
3984 printf("%f %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
3985 printf("%f %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
3986 printf("%f %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
3987 printf("%f %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
3988 printf("%f %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
3989 printf("%f %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
3990 printf("%f %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
3994 //------------------------------------------------------------------------
3995 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
3996 TString direction) {
3997 //-------------------------------------------------------------------
3998 // Propagate beyond beam pipe and correct for material
3999 // (material budget in different ways according to fUseTGeo value)
4000 // Add time if going outward (PropagateTo or PropagateToTGeo)
4001 //-------------------------------------------------------------------
4003 // Define budget mode:
4004 // 0: material from AliITSRecoParam (hard coded)
4005 // 1: material from TGeo in one step (on the fly)
4006 // 2: material from lut
4007 // 3: material from TGeo in one step (same for all hypotheses)
4020 if(fTrackingPhase.Contains("Clusters2Tracks"))
4021 { mode=3; } else { mode=1; }
4024 if(fTrackingPhase.Contains("Clusters2Tracks"))
4025 { mode=3; } else { mode=2; }
4031 if(fTrackingPhase.Contains("Default")) mode=0;
4033 Int_t index=fCurrentEsdTrack;
4035 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4036 Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4038 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4040 Double_t xOverX0,x0,lengthTimesMeanDensity;
4044 xOverX0 = AliITSRecoParam::GetdPipe();
4045 x0 = AliITSRecoParam::GetX0Be();
4046 lengthTimesMeanDensity = xOverX0*x0;
4047 lengthTimesMeanDensity *= dir;
4048 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4051 if (!t->PropagateToTGeo(xToGo,1)) return 0;
4054 if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4055 xOverX0 = fxOverX0Pipe;
4056 lengthTimesMeanDensity = fxTimesRhoPipe;
4057 lengthTimesMeanDensity *= dir;
4058 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4061 if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4062 if(fxOverX0PipeTrks[index]<0) {
4063 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4064 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4065 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4066 fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4067 fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4070 xOverX0 = fxOverX0PipeTrks[index];
4071 lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4072 lengthTimesMeanDensity *= dir;
4073 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4079 //------------------------------------------------------------------------
4080 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4082 TString direction) {
4083 //-------------------------------------------------------------------
4084 // Propagate beyond SPD or SDD shield and correct for material
4085 // (material budget in different ways according to fUseTGeo value)
4086 // Add time if going outward (PropagateTo or PropagateToTGeo)
4087 //-------------------------------------------------------------------
4089 // Define budget mode:
4090 // 0: material from AliITSRecoParam (hard coded)
4091 // 1: material from TGeo in steps of X cm (on the fly)
4092 // X = AliITSRecoParam::GetStepSizeTGeo()
4093 // 2: material from lut
4094 // 3: material from TGeo in one step (same for all hypotheses)
4107 if(fTrackingPhase.Contains("Clusters2Tracks"))
4108 { mode=3; } else { mode=1; }
4111 if(fTrackingPhase.Contains("Clusters2Tracks"))
4112 { mode=3; } else { mode=2; }
4118 if(fTrackingPhase.Contains("Default")) mode=0;
4120 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4122 Int_t shieldindex=0;
4123 if (shield.Contains("SDD")) { // SDDouter
4124 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4126 } else if (shield.Contains("SPD")) { // SPDouter
4127 rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4130 Error("CorrectForShieldMaterial"," Wrong shield name\n");
4134 // do nothing if we are already beyond the shield
4135 Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4136 if(dir<0 && rTrack > rToGo) return 1; // going outward
4137 if(dir>0 && rTrack < rToGo) return 1; // going inward
4141 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4143 Int_t index=2*fCurrentEsdTrack+shieldindex;
4145 Double_t xOverX0,x0,lengthTimesMeanDensity;
4150 xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4151 x0 = AliITSRecoParam::GetX0shield(shieldindex);
4152 lengthTimesMeanDensity = xOverX0*x0;
4153 lengthTimesMeanDensity *= dir;
4154 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4157 nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4158 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4161 if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4162 xOverX0 = fxOverX0Shield[shieldindex];
4163 lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4164 lengthTimesMeanDensity *= dir;
4165 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4168 if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4169 if(fxOverX0ShieldTrks[index]<0) {
4170 if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4171 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4172 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4173 fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4174 fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4177 xOverX0 = fxOverX0ShieldTrks[index];
4178 lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4179 lengthTimesMeanDensity *= dir;
4180 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4186 //------------------------------------------------------------------------
4187 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4189 Double_t oldGlobXYZ[3],
4190 TString direction) {
4191 //-------------------------------------------------------------------
4192 // Propagate beyond layer and correct for material
4193 // (material budget in different ways according to fUseTGeo value)
4194 // Add time if going outward (PropagateTo or PropagateToTGeo)
4195 //-------------------------------------------------------------------
4197 // Define budget mode:
4198 // 0: material from AliITSRecoParam (hard coded)
4199 // 1: material from TGeo in stepsof X cm (on the fly)
4200 // X = AliITSRecoParam::GetStepSizeTGeo()
4201 // 2: material from lut
4202 // 3: material from TGeo in one step (same for all hypotheses)
4215 if(fTrackingPhase.Contains("Clusters2Tracks"))
4216 { mode=3; } else { mode=1; }
4219 if(fTrackingPhase.Contains("Clusters2Tracks"))
4220 { mode=3; } else { mode=2; }
4226 if(fTrackingPhase.Contains("Default")) mode=0;
4228 Float_t dir = (direction.Contains("inward") ? 1. : -1.);
4230 Double_t r=fgLayers[layerindex].GetR();
4231 Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4233 Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4235 if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4237 Int_t index=6*fCurrentEsdTrack+layerindex;
4240 Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4243 // back before material (no correction)
4245 rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4246 if (!t->GetLocalXat(rOld,xOld)) return 0;
4247 if (!t->Propagate(xOld)) return 0;
4251 xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4252 lengthTimesMeanDensity = xOverX0*x0;
4253 lengthTimesMeanDensity *= dir;
4254 // Bring the track beyond the material
4255 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4258 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4259 if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4262 if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4263 xOverX0 = fxOverX0Layer[layerindex];
4264 lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4265 lengthTimesMeanDensity *= dir;
4266 // Bring the track beyond the material
4267 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4270 if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4271 if(fxOverX0LayerTrks[index]<0) {
4272 nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4273 if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4274 Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4275 ((1.-t->GetSnp())*(1.+t->GetSnp())));
4276 fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4277 fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4280 xOverX0 = fxOverX0LayerTrks[index];
4281 lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4282 lengthTimesMeanDensity *= dir;
4283 if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4290 //------------------------------------------------------------------------
4291 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4292 //-----------------------------------------------------------------
4293 // Initialize LUT for storing material for each prolonged track
4294 //-----------------------------------------------------------------
4295 fxOverX0PipeTrks = new Float_t[ntracks];
4296 fxTimesRhoPipeTrks = new Float_t[ntracks];
4297 fxOverX0ShieldTrks = new Float_t[ntracks*2];
4298 fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4299 fxOverX0LayerTrks = new Float_t[ntracks*6];
4300 fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4302 for(Int_t i=0; i<ntracks; i++) {
4303 fxOverX0PipeTrks[i] = -1.;
4304 fxTimesRhoPipeTrks[i] = -1.;
4306 for(Int_t j=0; j<ntracks*2; j++) {
4307 fxOverX0ShieldTrks[j] = -1.;
4308 fxTimesRhoShieldTrks[j] = -1.;
4310 for(Int_t k=0; k<ntracks*6; k++) {
4311 fxOverX0LayerTrks[k] = -1.;
4312 fxTimesRhoLayerTrks[k] = -1.;
4319 //------------------------------------------------------------------------
4320 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4321 //-----------------------------------------------------------------
4322 // Delete LUT for storing material for each prolonged track
4323 //-----------------------------------------------------------------
4324 if(fxOverX0PipeTrks) {
4325 delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4327 if(fxOverX0ShieldTrks) {
4328 delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4331 if(fxOverX0LayerTrks) {
4332 delete [] fxOverX0LayerTrks; fxOverX0LayerTrks = 0;
4334 if(fxTimesRhoPipeTrks) {
4335 delete [] fxTimesRhoPipeTrks; fxTimesRhoPipeTrks = 0;
4337 if(fxTimesRhoShieldTrks) {
4338 delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4340 if(fxTimesRhoLayerTrks) {
4341 delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4345 //------------------------------------------------------------------------
4346 void AliITStrackerMI::SetForceSkippingOfLayer() {
4347 //-----------------------------------------------------------------
4348 // Check if we are forced to skip layers
4349 // either we set to skip them in RecoParam
4350 // or they were off during data-taking
4351 //-----------------------------------------------------------------
4353 const AliEventInfo *eventInfo = GetEventInfo();
4355 for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4356 fForceSkippingOfLayer[l] = 0;
4358 if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4362 AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4363 AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4365 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4366 } else if(l==2 || l==3) {
4367 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4369 if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4375 //------------------------------------------------------------------------
4376 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4377 Int_t ilayer,Int_t idet) const {
4378 //-----------------------------------------------------------------
4379 // This method is used to decide whether to allow a prolongation
4380 // without clusters, because we want to skip the layer.
4381 // In this case the return value is > 0:
4382 // return 1: the user requested to skip a layer
4383 // return 2: track outside z acceptance
4384 //-----------------------------------------------------------------
4386 if (ForceSkippingOfLayer(ilayer)) return 1;
4388 Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4390 if (idet<0 && // out in z
4391 ilayer>innerLayCanSkip &&
4392 AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4393 // check if track will cross SPD outer layer
4394 Double_t phiAtSPD2,zAtSPD2;
4395 if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4396 if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4398 return 2; // always allow skipping, changed on 05.11.2009
4403 //------------------------------------------------------------------------
4404 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4405 Int_t ilayer,Int_t idet,
4406 Double_t dz,Double_t dy,
4407 Bool_t noClusters) const {
4408 //-----------------------------------------------------------------
4409 // This method is used to decide whether to allow a prolongation
4410 // without clusters, because there is a dead zone in the road.
4411 // In this case the return value is > 0:
4412 // return 1: dead zone at z=0,+-7cm in SPD
4413 // This method assumes that fSPDdetzcentre is ordered from -z to +z
4414 // return 2: all road is "bad" (dead or noisy) from the OCDB
4415 // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4416 // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4417 //-----------------------------------------------------------------
4419 // check dead zones at z=0,+-7cm in the SPD
4420 if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4421 Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4422 fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4423 fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4424 Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4425 fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4426 fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4427 for (Int_t i=0; i<3; i++)
4428 if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4429 AliDebug(2,Form("crack SPD %d track z %f %f %f %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4430 if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4434 // check bad zones from OCDB
4435 if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4437 if (idet<0) return 0;
4439 AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4442 Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4443 if (ilayer==0 || ilayer==1) { // ---------- SPD
4445 } else if (ilayer==2 || ilayer==3) { // ---------- SDD
4447 detSizeFactorX *= 2.;
4448 } else if (ilayer==4 || ilayer==5) { // ---------- SSD
4451 AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4452 if (detType==2) segm->SetLayer(ilayer+1);
4453 Float_t detSizeX = detSizeFactorX*segm->Dx();
4454 Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4456 // check if the road overlaps with bad chips
4458 if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4459 Float_t zlocmin = zloc-dz;
4460 Float_t zlocmax = zloc+dz;
4461 Float_t xlocmin = xloc-dy;
4462 Float_t xlocmax = xloc+dy;
4463 Int_t chipsInRoad[100];
4465 // check if road goes out of detector
4466 Bool_t touchNeighbourDet=kFALSE;
4467 if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4468 if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4469 if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4470 if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4471 AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
4473 // check if this detector is bad
4475 AliDebug(2,Form("lay %d bad detector %d",ilayer,idet));
4476 if(!touchNeighbourDet) {
4477 return 2; // all detectors in road are bad
4479 return 3; // at least one is bad
4483 Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4484 AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4485 if (!nChipsInRoad) return 0;
4487 Bool_t anyBad=kFALSE,anyGood=kFALSE;
4488 for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4489 if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4490 AliDebug(2,Form(" chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4491 if (det.IsChipBad(chipsInRoad[iCh])) {
4499 if(!touchNeighbourDet) {
4500 AliDebug(2,"all bad in road");
4501 return 2; // all chips in road are bad
4503 return 3; // at least a bad chip in road
4508 AliDebug(2,"at least a bad in road");
4509 return 3; // at least a bad chip in road
4513 if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4514 || !noClusters) return 0;
4516 // There are no clusters in road: check if there is at least
4517 // a bad SPD pixel or SDD anode or SSD strips on both sides
4519 Int_t idetInITS=idet;
4520 for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4522 if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4523 AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4526 //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4530 //------------------------------------------------------------------------
4531 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4532 const AliITStrackMI *track,
4533 Float_t &xloc,Float_t &zloc) const {
4534 //-----------------------------------------------------------------
4535 // Gives position of track in local module ref. frame
4536 //-----------------------------------------------------------------
4541 if(idet<0) return kTRUE; // track out of z acceptance of layer
4543 Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4545 Int_t lad = Int_t(idet/ndet) + 1;
4547 Int_t det = idet - (lad-1)*ndet + 1;
4549 Double_t xyzGlob[3],xyzLoc[3];
4551 AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4552 // take into account the misalignment: xyz at real detector plane
4553 if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4555 if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4557 xloc = (Float_t)xyzLoc[0];
4558 zloc = (Float_t)xyzLoc[2];
4562 //------------------------------------------------------------------------
4563 //------------------------------------------------------------------------
4564 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4566 // Method to be optimized further:
4567 // Aim: decide whether a track can be used for PlaneEff evaluation
4568 // the decision is taken based on the track quality at the layer under study
4569 // no information on the clusters on this layer has to be used
4570 // The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4571 // the cut is done on number of sigmas from the boundaries
4573 // Input: Actual track, layer [0,5] under study
4575 // Return: kTRUE if this is a good track
4577 // it will apply a pre-selection to obtain good quality tracks.
4578 // Here also you will have the possibility to put a control on the
4579 // impact point of the track on the basic block, in order to exclude border regions
4580 // this will be done by calling a proper method of the AliITSPlaneEff class.
4582 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4583 // return: Bool_t -> kTRUE if usable track, kFALSE if not usable.
4585 Int_t index[AliITSgeomTGeo::kNLayers];
4587 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4589 for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4590 index[k]=clusters[k];
4594 {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4595 AliITSlayer &layer=fgLayers[ilayer];
4596 Double_t r=layer.GetR();
4597 AliITStrackMI tmp(*track);
4599 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4600 Int_t ncl_out=0; Int_t ncl_in=0;
4601 for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4602 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4603 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4604 // if (tmp.GetClIndex(lay)>=0) ncl_out++;
4605 if(index[lay]>=0)ncl_out++;
4607 for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4608 AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4609 tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4610 if (index[lay]>=0) ncl_in++;
4612 Int_t ncl=ncl_out+ncl_in;
4613 Bool_t nextout = kFALSE;
4614 if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4615 else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4616 Bool_t nextin = kFALSE;
4617 if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4618 else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4619 // maximum number of missing clusters allowed in outermost layers
4620 if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4622 // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4623 if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4625 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout) return kFALSE;
4626 if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin) return kFALSE;
4627 if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4628 // if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff() && !tmp.GetConstrain()) return kFALSE;
4632 if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4633 Int_t idet=layer.FindDetectorIndex(phi,z);
4634 if(idet<0) { AliInfo(Form("cannot find detector"));
4637 // here check if it has good Chi Square.
4639 //propagate to the intersection with the detector plane
4640 const AliITSdetector &det=layer.GetDetector(idet);
4641 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4645 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4646 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4647 if(key>fPlaneEff->Nblock()) return kFALSE;
4648 Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4649 if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4651 // DEFINITION OF SEARCH ROAD FOR accepting a track
4653 Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4654 Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4655 Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4656 Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4657 // done for RecPoints
4659 // exclude tracks at boundary between detectors
4660 //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4661 Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4662 AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4663 AliDebug(2,Form("Local: track impact x=%f, z=%f",locx,locz));
4664 AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4665 if ( (locx-dx < blockXmn+boundaryWidth) ||
4666 (locx+dx > blockXmx-boundaryWidth) ||
4667 (locz-dz < blockZmn+boundaryWidth) ||
4668 (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4671 //------------------------------------------------------------------------
4672 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4674 // This Method has to be optimized! For the time-being it uses the same criteria
4675 // as those used in the search of extra clusters for overlapping modules.
4677 // Method Purpose: estabilish whether a track has produced a recpoint or not
4678 // in the layer under study (For Plane efficiency)
4680 // inputs: AliITStrackMI* track (pointer to a usable track)
4682 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4683 // traversed by this very track. In details:
4684 // - if a cluster can be associated to the track then call
4685 // AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4686 // - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4689 {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4690 AliITSlayer &layer=fgLayers[ilayer];
4691 Double_t r=layer.GetR();
4692 AliITStrackMI tmp(*track);
4696 if (!tmp.GetPhiZat(r,phi,z)) return;
4697 Int_t idet=layer.FindDetectorIndex(phi,z);
4699 if(idet<0) { AliInfo(Form("cannot find detector"));
4703 //propagate to the intersection with the detector plane
4704 const AliITSdetector &det=layer.GetDetector(idet);
4705 if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4709 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4711 Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4712 TMath::Sqrt(tmp.GetSigmaZ2() +
4713 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4714 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4715 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4716 Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4717 TMath::Sqrt(tmp.GetSigmaY2() +
4718 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4719 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4720 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4722 // road in global (rphi,z) [i.e. in tracking ref. system]
4723 Double_t zmin = tmp.GetZ() - dz;
4724 Double_t zmax = tmp.GetZ() + dz;
4725 Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4726 Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4728 // select clusters in road
4729 layer.SelectClusters(zmin,zmax,ymin,ymax);
4731 // Define criteria for track-cluster association
4732 Double_t msz = tmp.GetSigmaZ2() +
4733 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4734 AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4735 AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4736 Double_t msy = tmp.GetSigmaY2() +
4737 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4738 AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4739 AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4740 if (tmp.GetConstrain()) {
4741 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4742 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4744 msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4745 msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4747 msz = 1./msz; // 1/RoadZ^2
4748 msy = 1./msy; // 1/RoadY^2
4751 const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4753 Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4754 //Double_t tolerance=0.2;
4755 /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4756 idetc = cl->GetDetectorIndex();
4757 if(idet!=idetc) continue;
4758 //Int_t ilay = cl->GetLayer();
4760 if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4761 if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4763 Double_t chi2=tmp.GetPredictedChi2(cl);
4764 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4768 if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4770 AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4771 UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4772 if(key>fPlaneEff->Nblock()) return;
4773 Bool_t found=kFALSE;
4776 while ((cl=layer.GetNextCluster(clidx))!=0) {
4777 idetc = cl->GetDetectorIndex();
4778 if(idet!=idetc) continue;
4779 // here real control to see whether the cluster can be associated to the track.
4780 // cluster not associated to track
4781 if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4782 (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy > 1. ) continue;
4783 // calculate track-clusters chi2
4784 chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4785 // in particular, the error associated to the cluster
4786 //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4788 if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4790 if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4791 // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4792 // track->SetExtraModule(ilayer,idetExtra);
4794 if(!fPlaneEff->UpDatePlaneEff(found,key))
4795 AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4796 if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4797 Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values
4798 Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4799 Int_t cltype[2]={-999,-999};
4802 Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
4806 tr[2]=TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system
4807 tr[3]=TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is
4810 clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4811 clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4812 cltype[0]=layer.GetCluster(ci)->GetNy();
4813 cltype[1]=layer.GetCluster(ci)->GetNz();
4815 // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4816 // X->50/sqrt(12)=14 micron Z->450/sqrt(12)= 120 micron)
4817 // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4818 // It is computed properly by calling the method
4819 // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4821 //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4822 //if (tmp.PropagateTo(x,0.,0.)) {
4823 chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4824 AliCluster c(*layer.GetCluster(ci));
4825 c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4826 c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4827 //if (layer.GetCluster(ci)->GetGlobalCov(cov)) // by using this, instead, you got nominal cluster errors
4828 clu[2]=TMath::Sqrt(c.GetSigmaY2());
4829 clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4832 // Compute the angles between the track and the module
4833 // compute the angle "in phi direction", i.e. the angle in the transverse plane
4834 // between the normal to the module and the projection (in the transverse plane) of the
4836 // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
4837 Float_t tgl = tmp.GetTgl();
4838 Float_t phitr = tmp.GetSnp();
4839 phitr = TMath::ASin(phitr);
4840 Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
4842 Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
4843 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
4845 alpha = tmp.GetAlpha();
4846 Double_t phiglob = alpha+phitr;
4848 p[0] = TMath::Cos(phiglob);
4849 p[1] = TMath::Sin(phiglob);
4851 TVector3 pvec(p[0],p[1],p[2]);
4852 TVector3 normvec(rot[1],rot[4],rot[7]);
4853 Double_t angle = pvec.Angle(normvec);
4855 if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
4856 angle *= 180./TMath::Pi();
4859 TVector3 pt(p[0],p[1],0);
4860 TVector3 normt(rot[1],rot[4],0);
4861 Double_t anglet = pt.Angle(normt);
4863 Double_t phiPt = TMath::ATan2(p[1],p[0]);
4864 if(phiPt<0)phiPt+=2.*TMath::Pi();
4865 Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
4866 if(phiNorm<0) phiNorm+=2.*TMath::Pi();
4867 if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
4868 if(phiNorm>phiPt) anglet*=-1.;// pt-->normt clockwise: anglet>0
4869 if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
4870 anglet *= 180./TMath::Pi();
4872 AngleModTrack[2]=(Float_t) angle;
4873 AngleModTrack[0]=(Float_t) anglet;
4874 // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
4875 AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
4876 AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
4877 AngleModTrack[1]*=180./TMath::Pi(); // in degree
4879 fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);