X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITStrackerMI.cxx;h=59ddfa12f2a4bfad7ae79c9dda5fddaf71fd7acd;hb=70cb7fe43db2f2781930304a1dbb84513cbc42f0;hp=c6e1f169314605b6f5a6a840260d82889182e871;hpb=4edc3c8c04e812a7eeb19b4c1d17160f8ef39b3b;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITStrackerMI.cxx b/ITS/AliITStrackerMI.cxx index c6e1f169314..59ddfa12f2a 100644 --- a/ITS/AliITStrackerMI.cxx +++ b/ITS/AliITStrackerMI.cxx @@ -33,6 +33,7 @@ #include #include #include +#include #include "AliLog.h" #include "AliGeomManager.h" @@ -94,6 +95,7 @@ fDebugStreamer(0), fITSChannelStatus(0), fkDetTypeRec(0), fPlaneEff(0), +fSPDChipIntPlaneEff(0), fITSPid(0) { @@ -144,6 +146,7 @@ fDebugStreamer(0), fITSChannelStatus(0), fkDetTypeRec(0), fPlaneEff(0), +fSPDChipIntPlaneEff(0), fITSPid(0) { //-------------------------------------------------------------------- //This is the AliITStrackerMI constructor @@ -268,7 +271,11 @@ fITSPid(0) { Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff(); if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1) AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane)); - if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD(); + if (iplane<2) { + fPlaneEff = new AliITSPlaneEffSPD(); + fSPDChipIntPlaneEff = new Bool_t[AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip]; + for (UInt_t i=0; iGetReadPlaneEffFromOCDB()) @@ -354,6 +361,7 @@ AliITStrackerMI::~AliITStrackerMI() if(fITSChannelStatus) delete fITSChannelStatus; if(fPlaneEff) delete fPlaneEff; if(fITSPid) delete fITSPid; + if (fSPDChipIntPlaneEff) delete [] fSPDChipIntPlaneEff; } //------------------------------------------------------------------------ @@ -699,54 +707,58 @@ Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) { fTrackingPhase="PropagateBack"; Int_t nentr=event->GetNumberOfTracks(); // Info("PropagateBack", "Number of ESD tracks: %d\n", nentr); - + double bz0 = GetBz(); + const double kWatchStep=10.; // for larger steps watch arc vs segment difference + // Int_t ntrk=0; for (Int_t i=0; iGetTrack(i); - if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue; - if (esd->GetStatus()&AliESDtrack::kITSout) continue; - - AliITStrackMI *t = new AliITStrackMI(*esd); - - t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.)); - - ResetTrackToFollow(*t); - - /* - // propagate to vertex [SR, GSI 17.02.2003] - // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov - if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) { - if (fTrackToFollow.PropagateToVertex(event->GetVertex())) - fTrackToFollow.StartTimeIntegral(); - // from vertex to outside pipe - CorrectForPipeMaterial(&fTrackToFollow,"outward"); - }*/ // Start time integral and add distance from current position to vertex + if (esd->GetStatus()&AliESDtrack::kITSout) continue; + AliITStrackMI t(*esd); Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()}; - fTrackToFollow.GetXYZ(xyzTrk); + t.GetXYZ(xyzTrk); Double_t dst2 = 0.; - for (Int_t icoord=0; icoord<3; icoord++) { - Double_t di = xyzTrk[icoord] - xyzVtx[icoord]; - dst2 += di*di; + { + double dxs = xyzTrk[0] - xyzVtx[0]; + double dys = xyzTrk[1] - xyzVtx[1]; + double dzs = xyzTrk[2] - xyzVtx[2]; + // RS: for large segment steps d use approximation of cicrular arc s by + // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2) + // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature) + // Hence s^2/d^2 = (1+1/6 p^2)^2 + dst2 = dxs*dxs + dys*dys; + if (dst2 > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor + double crv = TMath::Abs(esd->GetC(bz0)); + double fcarc = 1.+crv*crv*dst2/6.; + dst2 *= fcarc*fcarc; + } + dst2 += dzs*dzs; } - fTrackToFollow.StartTimeIntegral(); - fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2)); + t.StartTimeIntegral(); + t.AddTimeStep(TMath::Sqrt(dst2)); + // + // transfer the time integral to ESD track + esd->SetStatus(AliESDtrack::kTIME); + Double_t times[10];t.GetIntegratedTimes(times); esd->SetIntegratedTimes(times); + esd->SetIntegratedLength(t.GetIntegratedLength()); + // + if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue; + t.SetExpQ(TMath::Max(0.8*t.GetESDtrack()->GetTPCsignal(),30.)); + ResetTrackToFollow(t); + // fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters(); - if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) { - if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) { - delete t; - continue; - } - fTrackToFollow.SetLabel(t->GetLabel()); - //fTrackToFollow.CookdEdx(); - CookLabel(&fTrackToFollow,0.); //For comparison only - fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout); - //UseClusters(&fTrackToFollow); - ntrk++; + if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,&t)) { + if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) continue; + fTrackToFollow.SetLabel(t.GetLabel()); + //fTrackToFollow.CookdEdx(); + CookLabel(&fTrackToFollow,0.); //For comparison only + fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout); + //UseClusters(&fTrackToFollow); + ntrk++; } - delete t; } AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr)); @@ -772,6 +784,13 @@ Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) { Int_t nentr=event->GetNumberOfTracks(); // Info("RefitInward", "Number of ESD tracks: %d\n", nentr); + // only for PlaneEff and in case of SPD (for FO studies) + if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() && + AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 && + AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) { + for (UInt_t i=0; iGetTrack(i); @@ -3072,7 +3091,7 @@ Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID else{ // Int_t tracks2[24], cluster[24]; - for (Int_t i=0;i=fTrackHypothesys.GetEntriesFast()) - fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10)); // TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex); if (!array) { array = new TObjArray(10); - fTrackHypothesys.AddAt(array,esdindex); + fTrackHypothesys.AddAtAndExpand(array,esdindex); } array->AddLast(track); } @@ -4786,6 +4803,7 @@ Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track, } } + if(zlocmin>zlocmax)return 0; Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax); AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad)); if (!nChipsInRoad) return 0; @@ -4867,7 +4885,7 @@ Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet, } //------------------------------------------------------------------------ //------------------------------------------------------------------------ -Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const { +Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer){ // // Method to be optimized further: // Aim: decide whether a track can be used for PlaneEff evaluation @@ -4903,19 +4921,19 @@ Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t AliITStrackMI tmp(*track); // require a minimal number of cluster in other layers and eventually clusters in closest layers - Int_t ncl_out=0; Int_t ncl_in=0; + Int_t nclout=0; Int_t nclin=0; for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay, tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ; - // if (tmp.GetClIndex(lay)>=0) ncl_out++; -if(index[lay]>=0)ncl_out++; + // if (tmp.GetClIndex(lay)>=0) nclout++; +if(index[lay]>=0)nclout++; } for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers AliDebug(2,Form("trak=%d lay=%d ; index=%d ESD label= %d",tmp.GetLabel(),lay, tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ; - if (index[lay]>=0) ncl_in++; + if (index[lay]>=0) nclin++; } - Int_t ncl=ncl_out+ncl_in; + Int_t ncl=nclout+nclin; Bool_t nextout = kFALSE; if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE ); @@ -4923,7 +4941,7 @@ if(index[lay]>=0)ncl_out++; if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE ); // maximum number of missing clusters allowed in outermost layers - if(ncl_outGetMaxMissingClustersOutPlaneEff()) + if(ncloutGetMaxMissingClustersOutPlaneEff()) return kFALSE; // maximum number of missing clusters allowed (both in innermost and in outermost layers) if(nclGetMaxMissingClustersPlaneEff()) @@ -4958,9 +4976,20 @@ if(index[lay]>=0)ncl_out++; // Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff(); Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff(); - Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()); // those are precisions in the tracking reference system - Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); // Use it also for the module reference system, as it is - // done for RecPoints + Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXFromBoundaryPlaneEff(); + Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZFromBoundaryPlaneEff(); + Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx; + // those are precisions in the tracking reference system + Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz; + // Use it also for the module reference system, as it is done for RecPoints + + if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigFrmBndPlaneEff()){ + if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) dx -= nsigx*TMath::Sqrt(tmp.GetSigmaY2()); + else dx -= distx; + + if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) dz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); + else dz -= distz; + } // exclude tracks at boundary between detectors //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff(); @@ -4972,6 +5001,28 @@ if(index[lay]>=0)ncl_out++; (locx+dx > blockXmx-boundaryWidth) || (locz-dz < blockZmn+boundaryWidth) || (locz+dz > blockZmx-boundaryWidth) ) return kFALSE; + + if(ilayer==0){ + const AliESDEvent *myesd = ((AliESDtrack*)tmp.GetESDtrack())->GetESDEvent(); + //The beam pipe + if (CorrectForPipeMaterial(&tmp,"inward")) { + const AliESDVertex* vtx = myesd->GetVertex(); + if(!vtx) return kFALSE; + Double_t ddz[2],cov[3]; + Double_t maxD=3.; + if(!tmp.PropagateToDCA(vtx,AliTracker::GetBz(),maxD,ddz,cov)) return kFALSE; + if(TMath::Abs(ddz[0])>=AliITSReconstructor::GetRecoParam()->GetDCACutPlaneEff()) return kFALSE; + + Double_t covar[6]; vtx->GetCovMatrix(covar); + Double_t p[2]={tmp.GetParameter()[0]-ddz[0],tmp.GetParameter()[1]-ddz[1]}; + Double_t c[3]={covar[2],0.,covar[5]}; + Double_t chi2= ((AliESDtrack*)tmp.GetESDtrack())->GetPredictedChi2(p,c); + if (chi2>AliITSReconstructor::GetRecoParam()->GetVertexChi2CutPlaneEff()) return kFALSE; // Use this to cut on chi^2 + + } else return kFALSE; + } + + return kTRUE; } //------------------------------------------------------------------------ @@ -5050,9 +5101,32 @@ void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilay msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC(); msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC(); } + + if(AliITSReconstructor::GetRecoParam()->GetSwitchOffStdSearchClusPlaneEff()){ + Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXSearchClusterPlaneEff(); + Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZSearchClusterPlaneEff(); + Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXSearchClusterPlaneEff(); + Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZSearchClusterPlaneEff(); + msy = nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx; + msz = nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz; + + if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigSrhClusPlaneEff()){ + if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) msy -= nsigx*TMath::Sqrt(tmp.GetSigmaY2()); + else msy -= distx; + + if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) msz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2()); + else msz -= distz; + } + + msy *= msy; + msz *= msz; + + } + + if(msz==0 || msy==0){AliWarning("UseTrackForPlaneEff: null search frame"); return;} + msz = 1./msz; // 1/RoadZ^2 msy = 1./msy; // 1/RoadY^2 -// const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1; Int_t idetc=-1; @@ -5099,13 +5173,36 @@ void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilay } if(!fPlaneEff->UpDatePlaneEff(found,key)) AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key)); + +// this for FO efficiency studies (only for SPD) // + UInt_t keyFO=999999; + Bool_t foundFO=kFALSE; + if(ilayer<2){ //ONLY SPD layers for FastOr studies + TBits mapFO = fkDetTypeRec->GetFastOrFiredMap(); + Int_t phase = (fEsd->GetBunchCrossNumber())%4; + if(!fSPDChipIntPlaneEff[key]){ + AliITSPlaneEffSPD spd; + keyFO = spd.SwitchChipKeyNumbering(key); + if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE; + keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1); + if(keyFOUpDatePlaneEff(foundFO,keyFO)) + AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO)); + } + } + + + if(fPlaneEff->GetCreateHistos()&& AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) { Float_t tr[4]={99999.,99999.,9999.,9999.}; // initialize to high values Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails) Int_t cltype[2]={-999,-999}; // and the module -Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below) + Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below) tr[0]=locx; tr[1]=locz; @@ -5175,15 +5272,19 @@ Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.; anglet *= 180./TMath::Pi(); - AngleModTrack[2]=(Float_t) angle; - AngleModTrack[0]=(Float_t) anglet; + angleModTrack[2]=(Float_t) angle; + angleModTrack[0]=(Float_t) anglet; // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90) - AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.)); - AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2 - AngleModTrack[1]*=180./TMath::Pi(); // in degree + angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.)); + angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2 + angleModTrack[1]*=180./TMath::Pi(); // in degree - fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack); + fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack); + + // For FO efficiency studies of SPD + if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack); } + if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE; return; } @@ -5232,4 +5333,3 @@ Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *stor } return nfound; } -