+//________________________________________________________________________________________________________
+Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
+{
+ // crude estimate of helix parameters, w/o errors and Eloss.
+ // Fast Riemann fit: Comp.Phy.Comm.131 (2000) 95
+ //
+ // if charge is not imposed (extQ==0) then it will be determined from the collision type
+ //
+ static TArrayD arrU,arrV,arrW;
+ double *parrW,*parrU,*parrV;
+ Bool_t eloss = IsELossON();
+ //
+ int np = fPntLast - fPntFirst + 1;
+ if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
+ //
+ const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov();
+ //
+ if (fPntLast>arrU.GetSize()) {
+ arrU.Set(2*fPntLast);
+ arrV.Set(2*fPntLast);
+ arrW.Set(2*fPntLast);
+ }
+ parrU = arrU.GetArray();
+ parrV = arrV.GetArray();
+ parrW = arrW.GetArray();
+ //
+ double uav=0,vav=0,wav=0,muu=0,muv=0,muw=0,mvv=0,mvw=0,mww=0;
+ int minRId = fPntFirst;
+ //
+ // get points span
+ double xmn=1e9,xmx=-1e9, ymn=1e9,ymx=-1e9;
+ for (int i=fPntFirst;i<=fPntLast;i++) {
+ parrW[i] = x[i]*x[i]+y[i]*y[i];
+ if (parrW[i]<parrW[minRId]) minRId = i; // point closest to origin
+ if (xmn>x[i]) xmn = x[i];
+ if (xmx<x[i]) xmx = x[i];
+ if (ymn>y[i]) ymn = y[i];
+ if (ymx<y[i]) ymx = y[i];
+ }
+ int minRId1 = minRId>fPntFirst ? fPntFirst:fPntFirst+1;
+ for (int i=fPntFirst;i<=fPntLast;i++) if (parrW[i]<parrW[minRId1] && i!=minRId) minRId1 = i;
+ //
+ double xshift = (xmx+xmn)/2 + 10*(ymx-ymn); // shift origin to have uniform weights
+ double yshift = (ymx+ymn)/2 - 10*(xmx-xmn);
+ // printf("X: %+e %+e Y: %+e %+e | shift: %+e %+e\n",xmn,xmx,ymn,ymx,xshift,yshift);
+ //
+ for (int i=fPntFirst;i<=fPntLast;i++) {
+ double xs = x[i] - xshift;
+ double ys = y[i] - yshift;
+ double w = xs*xs + ys*ys;
+ double scl = 1./(1.+w);
+ int i0 = i-fPntFirst;
+ wav += parrW[i0] = w*scl;
+ uav += parrU[i0] = xs*scl;
+ vav += parrV[i0] = ys*scl;
+ }
+ uav /= np; vav /= np; wav /= np;
+ //
+ for (int i=fPntFirst;i<=fPntLast;i++) {
+ //
+ // point next to closest
+ int i0 = i-fPntFirst;
+ if (parrW[i0]<parrW[minRId1-fPntFirst] && i!=minRId) minRId1 = i;
+ double u = parrU[i0] - uav;
+ double v = parrV[i0] - vav;
+ double w = parrW[i0] - wav;
+ muu += u*u;
+ muv += u*v;
+ muw += u*w;
+ mvv += v*v;
+ mvw += v*w;
+ mww += w*w;
+ }
+ muu/=np; muv/=np; muw/=np; mvv/=np; mvw/=np; mww/=np;
+ //
+ // find eigenvalues:
+ double trace3 = (muu + mvv + mww)/3.;
+ double muut = muu-trace3;
+ double mvvt = mvv-trace3;
+ double mwwt = mww-trace3;
+ double q = (muut*(mvvt*mwwt-mvw*mvw) - muv*(muv*mwwt-mvw*muw) + muw*(muv*mvw-mvvt*muw))/2;
+ double p = (muut*muut+mvvt*mvvt+mwwt*mwwt+2*(muv*muv+muw*muw+mvw*mvw))/6;
+ double dfpp = p*p*p-q*q;
+ dfpp = dfpp>0 ? TMath::Sqrt(dfpp)/q : 0;
+ double ph = TMath::ATan( dfpp )/3.;
+ if (ph<0) ph += TMath::Pi()/3;
+ p = p>0 ? TMath::Sqrt(p) : 0;
+ const double kSqrt3 = 1.73205080;
+ double snp = TMath::Sin(ph);
+ double csp = TMath::Cos(ph);
+ // double eg1 = trace3 + 2*p*csp;
+ double eg2 = trace3 - p*(csp+kSqrt3*snp); // smallest one
+ // double eg3 = trace3 - p*(csp-kSqrt3*snp);
+ // eigenvector for min.eigenvalue
+ muut = muu-eg2;
+ mvvt = mvv-eg2;
+ mwwt = mww-eg2;
+ double n0 = muv*mvw-muw*mvvt;
+ double n1 = muv*muw-mvw*muut;
+ double n2 = muut*mvvt-muv*muv;
+ // normalize to largest one
+ double nrm = TMath::Abs(n0);
+ if (nrm<TMath::Abs(n1)) nrm = TMath::Abs(n1);
+ if (nrm<TMath::Abs(n2)) nrm = TMath::Abs(n2);
+ n0/=nrm; n1/=nrm; n2/=nrm;
+ //
+ double cpar = -(uav*n0 + vav*n1 + wav*n2);
+ double xc = -n0/(cpar+n2)/2 + xshift;
+ double yc = -n1/(cpar+n2)/2 + yshift;
+ double rad = TMath::Sqrt(n0*n0+n1*n1-4*cpar*(cpar+n2))/2./TMath::Abs(cpar+n2);
+ //
+ // printf("Rad: %+e xc: %+e yc: %+e | X0: %+e Y0: %+e | X1: %+e Y1: %+e\n",rad,xc,yc, x[minRId],y[minRId],x[minRId1],y[minRId1]);
+
+ // linear circle fit --------------------------------------------------- <<<
+ //
+ // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
+ int sqb;
+ if (extQ) {
+ SetCharge(extQ);
+ sqb = fBz<0 ? -GetCharge():GetCharge();
+ }
+ else {
+ // determine the charge from the collision type and field sign
+ // the negative Q*B will have positive Vc x dir product Z component
+ // with Vc={-xc,-yc} : vector from circle center to the origin
+ // and V0 - track direction vector (take {0,-1,1} for cosmics)
+ // If Bz is not provided, assume positive Bz
+ if ( IsTypeCosmics() ) sqb = xc>0 ? -1:1;
+ else {
+ // track direction vector as a - diference between the closest and the next to closest to origin points
+ double v0X = x[minRId1] - x[minRId];
+ double v0Y = y[minRId1] - y[minRId];
+ sqb = (yc*v0X - xc*v0Y)<0 ? -1:1;
+ }
+ SetCharge( fBz<0 ? -sqb : sqb);
+ }
+ //
+ Double_t phi = TMath::ATan2(yc,xc);
+ if (sqb<0) phi += TMath::Pi();
+ if (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
+ else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
+ fParams[kPhi0] = phi;
+ fParams[kR0] = sqb<0 ? -rad:rad;
+ fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0];
+ //
+ // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
+ //
+ // find z-offset and dip + the parameter t of closest approach to hits - >>>
+ //
+ UInt_t hitLrPos=0; // pattern of hit layers at pos
+ UInt_t hitLrNeg=0; // and negative t's
+
+ Double_t ss=0,st=0,sz=0,stt=0,szt=0;
+ for (int i=fPntFirst;i<=fPntLast;i++) {
+ //
+ Double_t ze2 = cov[i*6 + kZZ];
+ Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z
+ if (fParams[kR0]<0) t += TMath::Pi();
+ if (t > TMath::Pi()) t -= TMath::Pi()*2;
+ else if (t <-TMath::Pi()) t += TMath::Pi()*2;
+ if (ze2<fgkAlmostZero) ze2 = 1E-8;
+ ze2 = 1./ze2;
+ ss += ze2;
+ st += t*ze2;
+ stt+= t*t*ze2;
+ sz += z[i]*ze2;
+ szt+= z[i]*t*ze2;
+ //
+ fCurT[i] = t; // parameter of the closest approach to the point
+ // printf("%d %+e %+e %+e %+e\n",i,x[i],y[i],z[i],t);
+ if (eloss) {
+ double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]);
+ int lr;
+ for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break;
+ if (lr<kMaxLrITS) {
+ if (t>0) hitLrPos |= (1<<lr); // set bit of the layer
+ else hitLrNeg |= (1<<lr); // set bit of the layer
+ }
+ }
+ }
+ double det = ss*stt - st*st;
+ if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence
+ fParams[kDZ] = sz/ss;
+ fParams[kDip] = 0;
+ }
+ else {
+ fParams[kDZ] = (sz*stt-st*szt)/det;
+ fParams[kDip] = -(ss*szt-st*sz)/det/fParams[kR0];
+ }
+ //
+ // find z-offset and dip + the parameter t of closest approach to hits - <<<
+ //
+ // fill info needed to account for ELoss ------------------------------- >>>
+ if (eloss) {
+ fNElsPnt = fPntLast - fPntFirst + 1;
+ //
+ // to account for the energy loss in the passive volumes, calculate the relevant t-parameters
+ double* tcur = fCurT + fPntFirst;
+ double* ecur = fElsDR+ fPntFirst;
+ //
+ for (int ilp=3;ilp--;) {
+ int id = fgkPassivLrITS[ilp];
+ double tp = GetHelixParAtR( fgkRLayITS[ id ] );
+ if (tp<0) continue; // does not hit this radius
+ //
+ tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp;
+ ecur[fNElsPnt] = fgRhoLITS[ id ];
+ fNElsPnt++;
+ // printf("Passive on lr %d %+e\n",ilp,tcur[fNElsPnt-1]);
+ //
+ if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics
+ tcur[fNElsPnt] = -tcur[fNElsPnt-1];
+ ecur[fNElsPnt] = ecur[fNElsPnt-1];
+ fNElsPnt++;
+ //printf("Passive* on lr %d %+e\n",ilp,-tcur[fNElsPnt-1]);
+ }
+ //
+ }
+ // check if some active layers did not miss the hit, treat them as passive
+ for (int ilp=6;ilp--;) {
+ int id = fgkActiveLrITS[ilp];
+ double tp = GetHelixParAtR( fgkRLayITS[ id ] );
+ if (tp<0) continue; // does not hit this radius
+ //
+ if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) {
+ tcur[fNElsPnt] = -tp;
+ ecur[fNElsPnt] = fgRhoLITS[ id ];
+ fNElsPnt++;
+ //printf("Missed on lr %d %+e\n",ilp,-tp);
+ }
+ //
+ if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) {
+ tcur[fNElsPnt] = tp;
+ ecur[fNElsPnt] = fgRhoLITS[ id ];
+ fNElsPnt++;
+ //printf("Missed* on lr %d %e\n",ilp,tp);
+ }
+ }
+ //
+ TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE); // index e-loss points in increasing order
+ // find the position of smallest positive t-param
+ for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break;
+ //
+ Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
+ Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy
+ Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass); // in the point of closest approach to beam
+ Double_t normS[3];
+ //
+ // Positive t-params: along the track direction for negative track, against for positive
+ Double_t pcur = ptot, ecurr = etot;
+ for (int ip=fFirstPosT;ip<fNElsPnt;ip++) {
+ int tID = fElsId[ip];
+ Double_t t = fCurT[ tID ];
+ //
+ if (tID>fPntLast) { // this is not a hit layer but passive layer
+ double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
+ xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
+ normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
+ normS[1] = -TMath::Sin(php);
+ normS[2] = 0;
+ }
+ else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module
+ fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
+ }
+ //
+ // negaive t-params: against the track direction for negative track, along for positive
+ pcur = ptot;
+ ecurr = etot;
+ for (int ip=fFirstPosT;ip--;) {
+ int tID = fElsId[ip];
+ Double_t t = fCurT[ tID ];
+ //
+ if (tID>fPntLast) { // this is not a hit layer but passive layer
+ double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
+ xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
+ normS[0] = -TMath::Cos(php); // normal to the cylinder at intersection point
+ normS[1] = -TMath::Sin(php);
+ normS[2] = 0;
+ }
+ else GetNormal(normS,fkPoints->GetCov()+tID*6); // vector normal to hit module
+ //
+ fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
+ }
+ }
+ // fill info needed to account for ELoss ------------------------------- <<<
+ //
+ return kTRUE;
+}
+
+/*