+void AliTPCcalibUnlinearity::AlignOROC(AliTPCseed *track, Int_t isec){
+ //
+ // fit the tracklet in OROC sepatately in Quadrant
+ //
+ //
+ if (isec<36) return;
+ const Int_t kMinClusterF=35;
+ const Int_t kMinClusterQ=10;
+ //
+ const Int_t kdrow1 =8; // rows to skip at the end
+ const Int_t kdrow0 =2; // rows to skip at beginning
+ const Float_t kedgey=3;
+ const Float_t kMaxDist=0.5;
+ const Float_t kMaxCorrY=0.1;
+ const Float_t kPRFWidth = 0.6; //cut 2 sigma of PRF
+ //
+ //
+ AliTPCROC * roc = AliTPCROC::Instance();
+ AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
+ const Double_t kXmean=roc->GetPadRowRadii(isec,53);
+ //
+ // full track fit parameters
+ //
+ TLinearFitter fyf(2,"pol1");
+ TLinearFitter fzf(2,"pol1");
+ TVectorD pyf(2), peyf(2),pzf(2), pezf(2);
+ Int_t nf=0;
+ //
+ // make full fit as reference
+ //
+ for (Int_t iter=0; iter<2; iter++){
+ fyf.ClearPoints();
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetDetector()!=isec) continue;
+ if (c->GetRow()<kdrow0) continue;
+ Double_t dx = c->GetX()-kXmean;
+ Double_t x[2]={dx, dx*dx};
+ if (iter==0 &&c->GetType()<0) continue;
+ if (iter==1){
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ }
+ fyf.AddPoint(x,c->GetY(),0.1);
+ fzf.AddPoint(x,c->GetZ(),0.1);
+ }
+ nf = fyf.GetNpoints();
+ if (nf<kMinClusterF) return; // not enough points - skip
+ fyf.Eval();
+ fyf.GetParameters(pyf);
+ fyf.GetErrors(peyf);
+ fzf.Eval();
+ fzf.GetParameters(pzf);
+ fzf.GetErrors(pezf);
+ }
+ //
+ // Make Fitters and params for 3 fitters
+ //
+ TLinearFitter *fitters[4];
+ Int_t npoints[4];
+ TVectorD params[4];
+ TVectorD errors[4];
+ Double_t chi2C[4];
+ for (Int_t i=0;i<4;i++) {
+ fitters[i] = new TLinearFitter(2,"pol1");
+ npoints[i]=0;
+ params[i].ResizeTo(2);
+ errors[i].ResizeTo(2);
+ }
+ //
+ // Update fitters
+ //
+ for (Int_t irow=kdrow0;irow<159-kdrow1;irow++) {
+ AliTPCclusterMI *c=track->GetClusterPointer(irow);
+ if (!c) continue;
+ if (c->GetDetector()!=isec) continue;
+ if (c->GetRow()<kdrow0) continue;
+ Double_t dx = c->GetX()-kXmean;
+ Double_t x[2]={dx, dx*dx};
+ Double_t yfit = pyf[0]+pyf[1]*dx;
+ Double_t dedge = c->GetX()*TMath::Tan(TMath::Pi()/18.)-TMath::Abs(yfit);
+ if (TMath::Abs(c->GetY()-yfit)>kMaxDist) continue;
+ if (dedge<kedgey) continue;
+ Double_t corrtrY =
+ corr->RPhiCOGCorrection(isec,c->GetRow(), c->GetPad(),
+ c->GetY(),yfit, c->GetZ(), pyf[1], c->GetMax(),2.5);
+ if (TMath::Abs(corrtrY)>kMaxCorrY) continue;
+ if (dx<0){
+ if (yfit<-kPRFWidth) fitters[0]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fitters[1]->AddPoint(x,c->GetY(),0.1);
+ }
+ if (dx>0){
+ if (yfit<-kPRFWidth) fitters[2]->AddPoint(x,c->GetY(),0.1);
+ if (yfit>kPRFWidth) fitters[3]->AddPoint(x,c->GetY(),0.1);
+ }
+ }
+
+ for (Int_t i=0;i<4;i++) {
+ npoints[i] = fitters[i]->GetNpoints();
+ if (npoints[i]>=kMinClusterQ){
+ fitters[i]->Eval();
+ Double_t chi2Fac = TMath::Sqrt(fitters[i]->GetChisquare()/(fitters[i]->GetNpoints()-2));
+ chi2C[i]=chi2Fac;
+ fitters[i]->GetParameters(params[i]);
+ fitters[i]->GetErrors(errors[i]);
+ // renormalize errors
+ errors[i][0]*=chi2Fac;
+ errors[i][1]*=chi2Fac;
+ }
+ }
+ //
+ // Fill fitters
+ //
+ for (Int_t i0=0;i0<4;i0++){
+ for (Int_t i1=0;i1<4;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ Int_t index0=i0*4+i1;
+ Double_t xy[1];
+ Double_t xfi[1];
+ xy[0] = pyf[1];
+ xfi[0] = pyf[1];
+ //
+ Double_t dy = params[i1][0]-params[i0][0];
+ Double_t sy = TMath::Sqrt(errors[i1][0]*errors[i1][0]+errors[i0][0]*errors[i0][0]);
+ Double_t dphi = params[i1][1]-params[i0][1];
+ Double_t sphi = TMath::Sqrt(errors[i1][1]*errors[i1][1]+errors[i0][1]*errors[i0][1]);
+ //
+ Int_t sector = isec-36;
+ //
+ if (TMath::Abs(dy/(sy+0.1))>5.) continue; // 5 sigma cut
+ if (TMath::Abs(dphi/(sphi+0.004))>5.) continue; // 5 sigma cut
+ TLinearFitter * fitterY,*fitterPhi;
+ fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*sector+index0);
+ if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
+ fitterY = (TLinearFitter*) fFitterQuadrantY.At(16*36+index0);
+ if (TMath::Abs(dy/(sy+0.1))<5.) fitterY->AddPoint(xy,dy,sy);
+ //
+ fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*sector+index0);
+ if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
+ fitterPhi = (TLinearFitter*) fFitterQuadrantPhi.At(16*36+index0);
+ if (TMath::Abs(dphi/(sphi+0.001))<5.) fitterPhi->AddPoint(xfi,dphi,sphi);
+ }
+ }
+ //
+ // Dump debug information
+ //
+ if (fStreamLevel>0){
+ TTreeSRedirector *cstream = GetDebugStreamer();
+ Int_t sector = isec-36;
+ if (cstream){
+ for (Int_t i0=0;i0<4;i0++){
+ for (Int_t i1=i0;i1<4;i1++){
+ if (i0==i1) continue;
+ if(npoints[i0]<kMinClusterQ) continue;
+ if(npoints[i1]<kMinClusterQ) continue;
+ (*cstream)<<"fitQ"<<
+ "run="<<fRun<< // run number
+ "event="<<fEvent<< // event number
+ "time="<<fTime<< // time stamp of event
+ "trigger="<<fTrigger<< // trigger
+ "mag="<<fMagF<< // magnetic field
+ "sec="<<sector<< // current sector from 0 to 36
+ "isec="<<isec<< // current sector
+ // full fit
+ "nf="<<nf<< // total number of points
+ "pyf.="<<&pyf<< // full OROC fit y
+ "pzf.="<<&pzf<< // full OROC fit z
+ // quadrant fit
+ "i0="<<i0<< // quadrant number
+ "i1="<<i1<<
+ "n0="<<npoints[i0]<< // number of points
+ "n1="<<npoints[i1]<<
+ "p0.="<<¶ms[i0]<< // parameters
+ "p1.="<<¶ms[i1]<<
+ "e0.="<<&errors[i0]<< // errors
+ "e1.="<<&errors[i1]<<
+ "chi0="<<chi2C[i0]<< // chi2s
+ "chi1="<<chi2C[i1]<<
+
+ "\n";
+ }
+ }
+ }
+ }
+}
+
+
+
+
+