+Bool_t AliTPCcalibLaser::FitDriftV(Float_t minFraction){
+ //
+ // Fit corrections to the drift velocity - linear approximation in the z and global y
+ //The transfromatiom from the drift time to the z position done in AliTPCTracnsform class
+ //
+ // Source of outlyers :
+ // 0. Track in the saturation - postpeak
+ // 1. gating grid close the part of the signal for first bundle
+
+ // The robust fit is performed in 2 itterations /robust fraction controlled by kFraction/
+ // 1. Robust fit is used in the itteration number 0
+ // only fraction of laser uted
+ // 2. Only the tracks close to the fit used in the second itteration
+ /*
+ Formulas:
+
+ z = s* (z0 - vd*(t-t0))
+
+ s - side -1 and +1
+ t0 - time 0
+ vd - nominal drift velocity
+ zs - miscalibrated position
+
+ zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
+ vr - relative change of the drift velocity
+ dzt - vd*dt
+ dr = zz0-s*z
+ ..
+ ==>
+ zs ~ z - s*vr*(z0-s*z)+s*dzt
+ --------------------------------
+ 1. Correction function vr constant:
+
+
+ dz = zs-z = -s*vr *(z0-s*z)+s*dzt
+ dzs/dl = dz/dl +s*s*vr*dz/dl
+ d(dz/dl) = vr*dz/dl
+ */
+ const Int_t knLaser = 336; //n laser tracks
+ const Float_t kFraction[2] = {0.70,0.95}; // robust fit fraction
+
+ const Float_t kSaturCut = 0.05; // remove saturated lasers - cut on fraction of saturated
+ const Float_t kDistCut = 3.; // distance sigma cut - 3 sigma
+ const Float_t kDistCutAbs = 1.; // absolute cut 1 cm
+ const Float_t kMinClusters = 40.; // minimal amount of the clusters
+ const Float_t kMinSignal = 2.5; // minimal mean height of the signal
+ const Float_t kChi2Cut = 1.0; // chi2 cut to accept drift fit
+ //
+ static TLinearFitter fdriftA(3,"hyp2");
+ static TLinearFitter fdriftC(3,"hyp2");
+ static TLinearFitter fdriftAC(4,"hyp3");
+ TVectorD fitA(3),fitC(3),fitAC(4);
+
+ AliTPCcalibDB* calib=AliTPCcalibDB::Instance();
+ AliTPCParam * tpcparam = calib->GetParameters();
+ //
+ // reset old data
+ //
+ for (Int_t id=0; id<336; id++) fFitZ[id]=0;
+ if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
+ if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
+ for (Int_t i=0;i<5; i++){
+ (*fFitCside)[i]=0;
+ (*fFitAside)[i]=0;
+ }
+ //
+ //
+ Float_t chi2A = 10;
+ Float_t chi2C = 10;
+ Float_t chi2AC = 10;
+ Int_t npointsA=0;
+ Int_t npointsC=0;
+ Int_t npointsAC=0;
+ Int_t nbA[4]={0,0,0,0};
+ Int_t nbC[4]={0,0,0,0};
+ TVectorD vecZM(336); // measured z potion of laser
+ TVectorD vecA(336); // accepted laser
+ TVectorD vecZF(336); // fitted position
+ TVectorD vecDz(336); // deltaZ
+ TVectorD vecZS(336); // surveyed position of laser
+ // additional variable to cut
+ TVectorD vecdEdx(336); // dEdx
+ TVectorD vecSy(336); // shape y
+ TVectorD vecSz(336); // shape z
+ //
+ //
+ for (Int_t id=0; id<336; id++){
+ Int_t reject=0;
+ AliTPCLaserTrack *ltrp =
+ (AliTPCLaserTrack*)AliTPCLaserTrack::GetTracks()->UncheckedAt(id);
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCseed * seed = (AliTPCseed*)fTracksTPC.At(id);
+ vecZM(id)= (param==0) ? 0:param->GetZ();
+ vecZS(id)= ltrp->GetZ();
+ vecDz(id)= 0;
+ vecA(id)=1;
+ vecdEdx(id)=(seed)?seed->GetdEdx():0;
+ vecSy(id) =(seed)?seed->CookShape(1):0;
+ vecSz(id) =(seed)?seed->CookShape(2):0;
+ //
+ fFitZ[id]=0;
+ if (!fTracksEsdParam.At(id)) reject|=1;
+ if (!param) continue;
+ if (!AcceptLaser(id)) reject|=2;
+ if ( fClusterSatur[id]>kSaturCut) reject|=4;
+ if ( fClusterCounter[id]<kMinClusters) reject|=8;
+ vecA(id)=reject;
+ if (reject>0) continue;
+ if (ltrp->GetSide()==0){
+ npointsA++;
+ nbA[ltrp->GetBundle()]++;
+ }
+ if (ltrp->GetSide()==1){
+ npointsC++;
+ nbC[ltrp->GetBundle()]++;
+ }
+ }
+ //
+ // reject "bad events"
+ //
+ Bool_t isOK=kTRUE;
+ Int_t naA=0;
+ Int_t naC=0;
+ if (npointsA<minFraction*0.5*knLaser && npointsC<minFraction*0.5*knLaser)
+ isOK=kFALSE;
+ for (Int_t i=0;i<4;i++){
+ //count accepted for all layers
+ if (nbA[i]>minFraction*0.5*0.25*knLaser) naA++;
+ if (nbC[i]>minFraction*0.5*0.25*knLaser) naC++;
+ }
+ if (naA<3 &&naC<3) isOK=kFALSE;
+ if (isOK==kFALSE) return kFALSE;
+
+ //
+ //
+ //
+ for (Int_t iter=0; iter<2; iter++){
+ fdriftA.ClearPoints();
+ fdriftC.ClearPoints();
+ fdriftAC.ClearPoints();
+ npointsA=0; npointsC=0; npointsAC=0;
+ //
+ for (Int_t id=0; id<336; id++){
+ if (!fTracksEsdParam.At(id)) continue;
+ if (!AcceptLaser(id)) continue;
+ if ( fClusterSatur[id]>kSaturCut) continue;
+ if ( fClusterCounter[id]<kMinClusters) continue;
+ AliESDtrack *track = (AliESDtrack*)fTracksEsd.At(id);
+ if (track->GetTPCsignal()<kMinSignal) continue;
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ Float_t sz = (ltrp->GetSide()==0) ? TMath::Sqrt(chi2A): TMath::Sqrt(chi2C);
+ if (npointsAC>0) sz =TMath::Sqrt(chi2AC);
+ if (iter>0 && TMath::Abs(fFitZ[id])>sz*kDistCut) continue;
+ if (iter>0 && TMath::Abs(fFitZ[id])>kDistCutAbs) continue;
+
+// // drift distance
+// Double_t zlength = tpcparam->GetZLength(0);
+// Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+// Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+ //
+ Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
+ Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
+ Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
+
+
+ Double_t xxx[2] = {ldrift,lxyz[1]*ldrift/(zlength*250.)};
+ if (iter==0 &<rp->GetBundle()==0) continue;
+ // skip bundle 0 - can be wrong in the 0 iteration
+ if (ltrp->GetSide()==0){
+ fdriftA.AddPoint(xxx,mdrift,1);
+ }else{
+ fdriftC.AddPoint(xxx,mdrift,1);
+ }
+ Double_t xxx2[3] = {ltrp->GetSide(), ldrift,lxyz[1]*ldrift/(zlength*250.)};
+ fdriftAC.AddPoint(xxx2,mdrift,1);
+ }
+ //
+ if (fdriftA.GetNpoints()>minFraction*0.5*knLaser){
+ //
+ fdriftA.Eval();
+ npointsA= fdriftA.GetNpoints();
+ chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
+ fdriftA.EvalRobust(kFraction[iter]);
+ fdriftA.GetParameters(fitA);
+ if (chi2A<kChi2Cut ||(*fFitAside)[0]==0 ) {
+ if (fFitAside->GetNoElements()<5) fFitAside->ResizeTo(5);
+ (*fFitAside)[0] = fitA[0];
+ (*fFitAside)[1] = fitA[1];
+ (*fFitAside)[2] = fitA[2];
+ (*fFitAside)[3] = fdriftA.GetNpoints();
+ (*fFitAside)[4] = chi2A;
+ }
+ }
+ if (fdriftC.GetNpoints()>minFraction*0.5*knLaser){
+ fdriftC.Eval();
+ npointsC= fdriftC.GetNpoints();
+ chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
+ fdriftC.EvalRobust(kFraction[iter]);
+ fdriftC.GetParameters(fitC);
+ if (chi2C<kChi2Cut||(*fFitCside)[0]==0) {
+ if (fFitCside->GetNoElements()<5) fFitCside->ResizeTo(5);
+ (*fFitCside)[0] = fitC[0];
+ (*fFitCside)[1] = fitC[1];
+ (*fFitCside)[2] = fitC[2];
+ (*fFitCside)[3] = fdriftC.GetNpoints();
+ (*fFitCside)[4] = chi2C;
+ }
+ }
+
+ if (fdriftAC.GetNpoints()>minFraction*knLaser &&npointsA>0.5*minFraction*knLaser&&npointsC>0.5*minFraction*knLaser){
+ fdriftAC.Eval();
+ npointsAC= fdriftAC.GetNpoints();
+ chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
+ fdriftAC.EvalRobust(kFraction[iter]);
+ fdriftAC.GetParameters(fitAC);
+ if (chi2AC<kChi2Cut||(*fFitACside)[0]==0) (*fFitACside) = fitAC;
+ (*fFitACside)[0] = fitAC[0];
+ (*fFitACside)[1] = fitAC[1];
+ (*fFitACside)[2] = fitAC[2];
+ (*fFitACside)[3] = fdriftAC.GetNpoints();
+ (*fFitACside)[4] = chi2AC;
+ }
+
+ for (Int_t id=0; id<336; id++){
+ if (!fTracksEsdParam.At(id)) continue;
+ //
+ AliExternalTrackParam *param=(AliExternalTrackParam*)fTracksEsdParam.At(id);
+ AliTPCLaserTrack *ltrp = ( AliTPCLaserTrack*)fTracksMirror.At(id);
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t lxyz[3];
+ Double_t lpxyz[3];
+ param->GetXYZ(xyz);
+ param->GetPxPyPz(pxyz);
+ ltrp->GetXYZ(lxyz);
+ ltrp->GetPxPyPz(lpxyz);
+ //Double_t zlength = tpcparam->GetZLength(0);
+ //Double_t ldrift = zlength-TMath::Abs(lxyz[2]);
+ //Double_t mdrift = zlength-TMath::Abs(xyz[2]);
+ Double_t zlength = (ltrp->GetSide()==0)? tpcparam->GetZLength(36): tpcparam->GetZLength(71);
+ Double_t ldrift = (ltrp->GetSide()==0)?zlength-lxyz[2]:lxyz[2]+zlength;
+ Double_t mdrift = (ltrp->GetSide()==0)?zlength-xyz[2]:xyz[2]+zlength;
+
+
+ Float_t fz =0;
+ if (ltrp->GetSide()==0){
+ fz = (fitA)[0]+(fitA)[1]*ldrift+(fitA)[2]*lxyz[1]*ldrift/(zlength*250.);
+ }else{
+ fz = (fitC)[0]+(fitC)[1]*ldrift+(fitC)[2]*lxyz[1]*ldrift/(zlength*250.);
+ }
+ if (npointsAC>10){
+ //fz = (fitAC)[0]+(fitAC)[1]*ltrp->GetSide()+(fitAC)[2]*ldrift+(fitAC)[3]*lxyz[1]*ldrift/(zlength*250.);
+ }
+ fFitZ[id]=mdrift-fz;
+ vecZF[id]=fz;
+ vecDz[id]=mdrift-fz;
+ }
+ if (fStreamLevel>0){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ TTimeStamp tstamp(fTime);
+ Float_t valuePressure0 = AliTPCcalibDB::GetPressure(tstamp,fRun,0);
+ Float_t valuePressure1 = AliTPCcalibDB::GetPressure(tstamp,fRun,1);
+ Double_t ptrelative0 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,0);
+ Double_t ptrelative1 = AliTPCcalibDB::GetPTRelative(tstamp,fRun,1);
+ Double_t temp0 = AliTPCcalibDB::GetTemperature(tstamp,fRun,0);
+ Double_t temp1 = AliTPCcalibDB::GetTemperature(tstamp,fRun,1);
+ TVectorD vecGoofie(20);
+ AliDCSSensorArray* goofieArray = AliTPCcalibDB::Instance()->GetGoofieSensors(fRun);
+ if (goofieArray)
+ for (Int_t isensor=0; isensor<goofieArray->NumSensors();isensor++){
+ AliDCSSensor *gsensor = goofieArray->GetSensor(isensor);
+ if (gsensor) vecGoofie[isensor]=gsensor->GetValue(tstamp);
+ }
+
+ if (cstream){
+ (*cstream)<<"driftvN"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ // Environment values
+ "press0="<<valuePressure0<<
+ "press1="<<valuePressure1<<
+ "pt0="<<ptrelative0<<
+ "pt1="<<ptrelative1<<
+ "temp0="<<temp0<<
+ "temp1="<<temp1<<
+ "vecGoofie.="<<&vecGoofie<<
+ //
+ //
+ "vecZM.="<<&vecZM<< // measured z position
+ "vecZS.="<<&vecZS<< // surveyed z position
+ "vecZF.="<<&vecZF<< // fitted z position
+ "vecDz.="<<&vecDz<< // fitted z position
+ "vecA.="<<&vecA<< // accept laser flag
+ "vecdEdx.="<<&vecdEdx<< // dEdx - to cut on
+ "vecSy.="<<&vecSy<< // shape y - to cut on
+ "vecSz.="<<&vecSz<< // shape z - to cut on
+ //
+ "iter="<<iter<<
+ "driftA.="<<fFitAside<<
+ "driftC.="<<fFitCside<<
+ "driftAC.="<<fFitACside<<
+ "chi2A="<<chi2A<<
+ "chi2C="<<chi2C<<
+ "chi2AC="<<chi2AC<<
+ "nA="<<npointsA<<
+ "nC="<<npointsC<<
+ "nAC="<<npointsAC<<
+ "\n";
+ /*
+ //
+ variables to check in debug mode:
+ //
+ chainDriftN->SetAlias("driftS","250-abs(vecZS.fElements)");
+ chainDriftN->SetAlias("driftM","250-abs(vecZM.fElements)");
+ chainDriftN->SetAlias("driftF","vecZF.fElements");
+ chainDriftN->SetAlias("deltaZ","driftF-driftM"); //deltaZ
+ TCut cutA="iter==1&&sqrt(chi2A)<0.1&&vecZS.fElements>0&&vecA.fElements==0";
+ TCut cutC="iter==1&&sqrt(chi2C)<0.1&&vecZS.fElements<0&&vecA.fElements==0";
+
+ */
+ }
+ }
+ }
+ return kTRUE;
+}