X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=ITS%2FAliITSRiemannFit.cxx;h=6a08853e4cd26db521478eaf776b12c650bd2908;hb=182e22d97a879a1bad8db77feaec21544835e3d2;hp=5c9ea1defa5b2546e905ed426ed92f38582033ae;hpb=5d12ce3847b64498d417a697c24fb86f715f2752;p=u%2Fmrichter%2FAliRoot.git diff --git a/ITS/AliITSRiemannFit.cxx b/ITS/AliITSRiemannFit.cxx index 5c9ea1defa5..6a08853e4cd 100644 --- a/ITS/AliITSRiemannFit.cxx +++ b/ITS/AliITSRiemannFit.cxx @@ -13,24 +13,26 @@ * provided "as is" without express or implied warranty. * * * * * - * ///////////////////////////////////////////////////////////////////// * - * * - * This class performs a fast fit of helices going through the <=6 * - * points of the ITS, with the goal of studying tracking and * - * vertexing performances. * - * Generated kinematics is used to take into account different weights * - * associated to points in different layers (with different multiple * - * scattering-originated errors). * - * * - * Based on the work by A. Strandlie, R. Fruhwirth * - * * - * First implementation by N. Bustreo, R. Turrisi - July 2000 * - * * - * Further modifications by A. Dainese, R. Turrisi * - * * - * Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it * - * * - * **************************************************************************/ + **************************************************************************/ + +// +// * +// This class performs a fast fit of helices going through the <=6 * +// points of the ITS, with the goal of studying tracking and * +// vertexing performances. * +// Generated kinematics is used to take into account different weights * +// associated to points in different layers (with different multiple * +// scattering-originated errors). * +// * +// Based on the work by A. Strandlie, R. Fruhwirth * +// * +// First implementation by N. Bustreo, R. Turrisi - July 2000 * +// * +// Further modifications by A. Dainese, R. Turrisi * +// * +// Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it * +// * +// ************************************************************************ // // // Modified November, 7th 2001 by Rosario Turrisi @@ -49,28 +51,30 @@ // // "PROPER WEIGHTS": (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2) // + + + #include +#include "AliITS.h" #include "AliITSRiemannFit.h" #include "AliRun.h" #include "TClonesArray.h" #include "stdio.h" #include "stdlib.h" #include "Riostream.h" -#include "TH2.h" #include "TMath.h" #include "TF1.h" #include "TGraphErrors.h" -#include "TMinuit.h" -#include "TCanvas.h" #include "TStyle.h" -#include "TRandom.h" #include "TParticle.h" -#include "TFile.h" +#include "TTree.h" +#include "TVector3.h" #include "AliITSRecPoint.h" #include "AliITSgeom.h" #include "AliITSmodule.h" #include "AliMC.h" + ClassImp(AliITSRiemannFit) @@ -91,8 +95,23 @@ AliITSRiemannFit::AliITSRiemannFit() { for(Int_t i=0;i<6;i++)fPLay[i] = 0; } -//---------------------------------------------------------------------- +//______________________________________________________________________ +AliITSRiemannFit::AliITSRiemannFit(const AliITSRiemannFit &rf) : TObject(rf) { + // Copy constructor + // Copies are not allowed. The method is protected to avoid misuse. + Error("AliITSRiemannFit","Copy constructor not allowed\n"); +} + +//______________________________________________________________________ +AliITSRiemannFit& AliITSRiemannFit::operator=(const AliITSRiemannFit& /* rf */){ + // Assignment operator + // Assignment is not allowed. The method is protected to avoid misuse. + Error("= operator","Assignment operator not allowed\n"); + return *this; +} + +//______________________________________________________________________ AliITSRiemannFit::~AliITSRiemannFit() { /////////////////////////////////////////////////////////// // Default destructor. @@ -129,24 +148,50 @@ AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks) { // test erase // fspdi = 0; // fspdo = 0; - Point_tl *first = new Point_tl[fSizeEvent]; - Point_tl **PointRecs = new Point_tl*[fSizeEvent]; + AliPointtl *first = new AliPointtl[fSizeEvent]; + AliPointtl **pointRecs = new AliPointtl*[fSizeEvent]; for(Int_t i=0;i<6;i++)fPLay[i] = 0; for(Int_t j=0;jphi = phi; - Points[i]->eta = -0.5*tan(0.5*TMath::ATan2(r,z)); - Points[i]->fx = x; - Points[i]->fy = y; - Points[i]->fz = z; - Points[i]->fdx = error[0]; - Points[i]->fdy = error[1]; - Points[i]->fdz = error[2]; - Points[i]->fr = r; - Points[i]->track = track; - Points[i]->lay = id[0], - Points[i]->lad = id[1]; - Points[i]->det = id[2]; - Points[i]->fMomentum = PE; - Points[i]->fOrigin = OT; - Points[i]->fPt = sqrt(PE.X()*PE.X()+PE.Y()*PE.Y()); - Points[i]->fCode = code; - Points[i]->fName = name; - Points[i]->vertexPhi = phiorigin; + if(phi<0.0) phi += pPI2; + Points[i]->SetPhi(phi); + Points[i]->SetEta(-0.5*tan(0.5*TMath::ATan2(r,z))); + Points[i]->SetX(x); + Points[i]->SetY(y); + Points[i]->SetZ(z); + Points[i]->SetdX(error[0]); + Points[i]->SetdY(error[1]); + Points[i]->SetdZ(error[2]); + Points[i]->SetR(r); + Points[i]->SetTrack(track); + Points[i]->SetLay(id[0]); + Points[i]->SetLad(id[1]); + Points[i]->SetDet(id[2]); + Points[i]->SetMomentum(&pE); + Points[i]->SetOrigin(&oT); + Points[i]->SetPt(sqrt(pE.X()*pE.X()+pE.Y()*pE.Y())); + Points[i]->SetCode(code); + Points[i]->SetName(name); + Points[i]->SetVertexPhi(phiorigin); index++; return; @@ -192,8 +237,8 @@ void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS, TParticle *part; AliITSgeom *gm = (AliITSgeom*)ITS->GetITSgeom(); //get pointer to modules array - TObjArray *ITSmodules = ITS->GetModules(); - Int_t nmodules=ITSmodules->GetEntriesFast(); + TObjArray *iTSmodules = ITS->GetModules(); + Int_t nmodules=iTSmodules->GetEntriesFast(); printf("nmodules = %d \n",nmodules); // Get the points from points file AliITSmodule *itsModule; @@ -201,44 +246,35 @@ void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS, Stat_t nent; AliITSRecPoint *recp; nent=TR->GetEntries(); - TClonesArray *ITSrec = ITS->RecPoints(); + TClonesArray *iTSrec = ITS->RecPoints(); - Int_t TotRP=0; + Int_t totRP=0; for (mod=0; modAt(mod); + itsModule=(AliITSmodule*)iTSmodules->At(mod); ITS->ResetRecPoints(); TR->GetEvent(mod); - Int_t nrecp = ITSrec->GetEntries(); + Int_t nrecp = iTSrec->GetEntries(); if(!nrecp) continue; - TotRP += nrecp; + totRP += nrecp; } - Int_t iMAX = TotRP; + Int_t iMAX = totRP; fPrimaryTracks = ntracks; fParticles = nparticles; - Point_tl *global = new Point_tl[iMAX]; - fPointRecs = new Point_tl*[iMAX]; + AliITSRiemannFit::AliPointtl *global = new AliPointtl[iMAX]; + fPointRecs = new AliITSRiemannFit::AliPointtl*[iMAX]; // - // test erase -// Point_tl *first = new Point_tl[iMAX]; -// Point_tl *second = new Point_tl[iMAX]; -// fspdi = new Point_tl*[iMAX]; -// fspdo = new Point_tl*[iMAX]; for(Int_t j=0;jAt(mod); + itsModule=(AliITSmodule*)iTSmodules->At(mod); ITS->ResetRecPoints(); TR->GetEvent(mod); - Int_t nrecp = ITSrec->GetEntries(); + Int_t nrecp = iTSrec->GetEntries(); if (!nrecp) continue; itsModule->GetID(layer,ladder,detector); for (irec=0;irecUncheckedAt(irec); + recp = (AliITSRecPoint*)iTSrec->UncheckedAt(irec); track=recp->fTracks[0]; if(track <0 ) continue; xcluster=recp->GetX(); // x on cluster zcluster=recp->GetZ(); // z on cluster part = (TParticle*) gAlice->GetMCApp()->Particle(track); - part->ProductionVertex(OT); // set the vertex - part->Momentum(PE); // set the vertex momentum + part->ProductionVertex(oT); // set the vertex + part->Momentum(pE); // set the vertex momentum name = part->GetName(); + Char_t nam2[50]; + sprintf(nam2,"%s",name); code = part->GetPdgCode(); - Phi = part->Phi(); + pPhi = part->Phi(); id[0]=layer; id[1]=ladder; id[2]=detector; @@ -272,31 +310,31 @@ void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS, locals[0]=xcluster; // x on cluster locals[1]=0.0; // y on cluster locals[2]=zcluster; // z on cluster - locals_error[0]=sqrt(recp->GetSigmaX2()); - locals_error[1]=0.0; - locals_error[2]=sqrt(recp->GetSigmaZ2()); - locals_plus[0]=xcluster+sqrt(recp->GetSigmaX2()); // x on cluster - if(layer==1||layer==2) locals_plus[1]=0.0150/2; // y on cluster - else if(layer==3||layer==4) locals_plus[1]=0.0280/2; // y on cluster - else if(layer==5||layer==6) locals_plus[1]=0.0300/2; // y on cluster - locals_plus[2]=zcluster+sqrt(recp->GetSigmaZ2()); // z on cluster - locals_minus[0]=xcluster-sqrt(recp->GetSigmaX2()); // x on cluster - locals_minus[1]=0.0; // y on cluster - locals_minus[2]=zcluster-sqrt(recp->GetSigmaZ2()); // z on cluster + localserror[0]=sqrt(recp->GetSigmaX2()); + localserror[1]=0.0; + localserror[2]=sqrt(recp->GetSigmaZ2()); + localsplus[0]=xcluster+sqrt(recp->GetSigmaX2()); // x on cluster + if(layer==1||layer==2) localsplus[1]=0.0150/2; // y on cluster + else if(layer==3||layer==4) localsplus[1]=0.0280/2; // y on cluster + else if(layer==5||layer==6) localsplus[1]=0.0300/2; // y on cluster + localsplus[2]=zcluster+sqrt(recp->GetSigmaZ2()); // z on cluster + localsminus[0]=xcluster-sqrt(recp->GetSigmaX2()); // x on cluster + localsminus[1]=0.0; // y on cluster + localsminus[2]=zcluster-sqrt(recp->GetSigmaZ2()); // z on cluster gm->LtoG(layer,ladder,detector,locals,xpoint); - gm->LtoG(layer,ladder,detector,locals_plus,error_plus); - gm->LtoG(layer,ladder,detector,locals_minus,error_minus); - global_error[0]=0.5*TMath::Abs(error_plus[0]-error_minus[0]); - global_error[1]=0.5*TMath::Abs(error_plus[1]-error_minus[1]); - global_error[2]=0.5*TMath::Abs(error_plus[2]-error_minus[2]); + gm->LtoG(layer,ladder,detector,localsplus,errorPlus); + gm->LtoG(layer,ladder,detector,localsminus,errorMinus); + globalError[0]=0.5*TMath::Abs(errorPlus[0]-errorMinus[0]); + globalError[1]=0.5*TMath::Abs(errorPlus[1]-errorMinus[1]); + globalError[2]=0.5*TMath::Abs(errorPlus[2]-errorMinus[2]); if(trackEta())<=1.0) ieta++; if(TMath::Abs(part->Eta())<=0.5) ieta2++; } if(!(id[0]==idold[0]&&id[1]==idold[1]&& id[2]==idold[2]&&id[3]==idold[3])) { - FillPoints(fPointRecs,num,xpoint,global_error,PE,OT,id,track,name,code,Phi); + FillPoints(fPointRecs,num,xpoint,globalError,pE,oT,id,track,nam2,code,pPhi); // // test erase switch (idold[0]) { @@ -320,11 +358,11 @@ void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS, break; } // if(idold[0]==1){ -// FillPoints(fspdi,nspdi,xpoint,global_error,PE,OT,id,track,name,code,Phi); +// FillPoints(fspdi,nspdi,xpoint,globalError,pE,oT,id,track,name,code,pPhi); // } // if(idold[0]==2){ -// FillPoints(fspdo,nspdo,xpoint,global_error,PE,OT,id,track,name,code,Phi); +// FillPoints(fspdo,nspdo,xpoint,globalError,pE,oT,id,track,name,code,pPhi); // } // if(idold[0]==3){ // nsddi++; @@ -365,27 +403,27 @@ void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS, /////////////////////////////////////////////////////////// // Functions for sorting the fPointRecs array /////////////////////////////////////////////////////////// -Bool_t SortZ(const Point_tl *s1,const Point_tl *s2){ +Bool_t SortZ(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){ // Z sorting function for qsort. Float_t a; - a = s1->fz - s2->fz; + a = s1->GetZ() - s2->GetZ(); if(a<0.0) return kTRUE; if(a>0.0) return kFALSE; return kFALSE; } -Bool_t SortTrack(const Point_tl *s1,const Point_tl *s2){ +Bool_t SortTrack(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){ // track sorting function for qsort. Float_t a; - a = s1->track - s2->track; + a = s1->GetTrack() - s2->GetTrack(); if(a<0.0) return kTRUE; if(a>0.0) return kFALSE; return kFALSE; } -void hpsortTrack(Point_tl **ra,Int_t n){ +void hpsortTrack(AliITSRiemannFit::AliPointtl **ra,Int_t n){ Int_t i,ir,j,l; - Point_tl *rra; + AliITSRiemannFit::AliPointtl *rra; if(n<2) return; @@ -417,9 +455,9 @@ void hpsortTrack(Point_tl **ra,Int_t n){ ra[i] = rra; } // end for ever } -void hpsortZ(Point_tl **ra,Int_t n){ +void hpsortZ(AliITSRiemannFit::AliPointtl **ra,Int_t n){ Int_t i,ir,j,l; - Point_tl *rra; + AliITSRiemannFit::AliPointtl *rra; if(n<2) return; @@ -488,9 +526,9 @@ void AliITSRiemannFit::WritePoints(void) { ///////////////////////////////////////////////////////////////////// FILE *ascii= fopen("AsciiPoints.dat","w"); for(Int_t i=0;ilay, - fPointRecs[i]->track,fPointRecs[i]->fx,fPointRecs[i]->fy, - fPointRecs[i]->fz); + fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->GetLay(), + fPointRecs[i]->GetTrack(),fPointRecs[i]->GetX(), + fPointRecs[i]->GetY(),fPointRecs[i]->GetZ()); } fclose(ascii); return; @@ -504,10 +542,11 @@ void AliITSRiemannFit::ReadPoints(void) { hpsortTrack(fPointRecs,fPoints); for(Int_t i=0;ilay,fPointRecs[i]->track,fPointRecs[i]->fx, - fPointRecs[i]->fy,fPointRecs[i]->fz,fPointRecs[i]->fOrigin.X(), - fPointRecs[i]->fOrigin.Y(),fPointRecs[i]->fOrigin.Z(), - fPointRecs[i]->fPt,fPointRecs[i]->fName); + i,fPointRecs[i]->GetLay(),fPointRecs[i]->GetTrack(), + fPointRecs[i]->GetX(),fPointRecs[i]->GetY(), + fPointRecs[i]->GetZ(),fPointRecs[i]->GetOrigin()->X(), + fPointRecs[i]->GetOrigin()->Y(),fPointRecs[i]->GetOrigin()->Z(), + fPointRecs[i]->GetPt(),fPointRecs[i]->GetName()); return; } //----------------------------------------------------------------------- @@ -521,14 +560,14 @@ Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c, /// returns x1 , x2 , x3 //////////////////////////////////////// - Double_t Q = ((a*a - 3*b)/9); - Double_t R = ((2*a*a*a - 9*a*b +27*c)/54); + Double_t qQ = ((a*a - 3*b)/9); + Double_t rR = ((2*a*a*a - 9*a*b +27*c)/54); Double_t theta; - Double_t F = -2*sqrt(Q); + Double_t aF = -2*sqrt(qQ); Double_t g = a/3; - Double_t PI2 = TMath::Pi()*2; + Double_t pPI2 = TMath::Pi()*2; - if( R*R>Q*Q*Q ) { + if( rR*rR>qQ*qQ*qQ ) { cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<X()*From[i]->X()+From[i]->Y()*From[i]->Y()); - Theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X()); - if(Theta[i]<0) Theta[i]+=PI2; - x = R[i]*cos(Theta[i])/(1+R[i]*R[i]); - y = R[i]*sin(Theta[i])/(1+R[i]*R[i]); - z = R[i]*R[i]/(1+R[i]*R[i]); + rR[i] = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y()); + theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X()); + if(theta[i]<0) theta[i]+=pPI2; + x = rR[i]*cos(theta[i])/(1+rR[i]*rR[i]); + y = rR[i]*sin(theta[i])/(1+rR[i]*rR[i]); + z = rR[i]*rR[i]/(1+rR[i]*rR[i]); To[i]->SetXYZ(x,y,z); } - delete[] R; - delete[] Theta; + delete[] rR; + delete[] theta; return; } @@ -575,7 +614,7 @@ void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) { Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t omega, Double_t &thu0, Double_t &thv0, Double_t &phi, TVector2 &zData, TVector3 &zError, - Double_t &CorrLin){ + Double_t &corrLin){ /////////////////////////////////////////////////////////////////////// // Fit the points in the (z,s) plane - helix 3rd equation // @@ -628,7 +667,7 @@ Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t ome if ( TMath::Abs(direction) != (npoints-1) ) {return 11;} TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez); - fitHist->Fit("pol1","Q"); + fitHist->Fit("pol1","qQ"); z0 = fitHist->GetFunction("pol1")->GetParameter(0); vpar = fitHist->GetFunction("pol1")->GetParameter(1); ez0 = fitHist->GetFunction("pol1")->GetParError(0); @@ -637,32 +676,32 @@ Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t ome zData.Set(z0,vpar); zError.SetXYZ(ez0,evpar,chisquare); - Double_t Sigmas=0.; - Double_t Sigmaz=0.; - Double_t Avs=0.; - Double_t Avz=0.; - Double_t Avsz=0.; + Double_t sigmas=0.; + Double_t sigmaz=0.; + Double_t avs=0.; + Double_t avz=0.; + Double_t avsz=0.; for(Int_t j = 0; j < npoints; j++) { - Avs += s[j]; - Avz += z[j]; - Avsz += s[j]*z[j]; + avs += s[j]; + avz += z[j]; + avsz += s[j]*z[j]; } - Avs /= (Double_t)npoints; - Avz /= (Double_t)npoints; - Avsz /= (Double_t)npoints; + avs /= (Double_t)npoints; + avz /= (Double_t)npoints; + avsz /= (Double_t)npoints; for(Int_t l = 0; l < npoints; l++) { - Sigmas += (s[l]-Avs)*(s[l]-Avs); - Sigmaz += (z[l]-Avz)*(z[l]-Avz); + sigmas += (s[l]-avs)*(s[l]-avs); + sigmaz += (z[l]-avz)*(z[l]-avz); } - Sigmas /=(Double_t)npoints; - Sigmaz /=(Double_t)npoints; + sigmas /=(Double_t)npoints; + sigmaz /=(Double_t)npoints; - Sigmas = sqrt(Sigmas); - Sigmaz = sqrt(Sigmaz); + sigmas = sqrt(sigmas); + sigmaz = sqrt(sigmaz); - CorrLin = (Avsz-Avs*Avz)/(Sigmas*Sigmaz); + corrLin = (avsz-avs*avz)/(sigmas*sigmaz); delete [] z; delete [] x; @@ -707,37 +746,37 @@ Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Doubl u0 = -9999.9999; // parameters of helix - strange value... v0 = -9999.9999; // parameters of helix - strange value... rho = -9999.9999; // parameters of helix -unphysical strange value... - Int_t Layer = 0; + Int_t pLayer = 0; const Char_t* name = 0; Int_t i=0,k=0; Int_t iMAX = 50; - Int_t N = 0; + Int_t nN = 0; Int_t npl[6]={0,0,0,0,0,0}; - Double_t P = sqrt(Px*Px+Py*Py+Pz*Pz); - Double_t Pt = sqrt(Px*Px+Py*Py); + Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz); + Double_t pt = sqrt(Px*Px+Py*Py); TVector3 zError; TVector2 zData; - Double_t CorrLin; + Double_t corrLin; TVector3 *ori = new TVector3[iMAX]; TVector3 **original = new TVector3*[iMAX]; TVector3 *rie = new TVector3[iMAX]; - TVector3 **Riemann = new TVector3*[iMAX]; + TVector3 **riemann = new TVector3*[iMAX]; TVector3 *err = new TVector3[iMAX]; TVector3 **errors = new TVector3*[iMAX]; TVector3 *linerr = new TVector3[iMAX]; TVector3 **linerrors = new TVector3*[iMAX]; - //PH Double_t Weight[iMAX]; - Double_t * Weight = new Double_t[iMAX]; + //PH Double_t weight[iMAX]; + Double_t * weight = new Double_t[iMAX]; for(i=0;iSetXYZ(9999,9999,9999); - Double_t A11,A12,A13,A21,A22,A23,A31,A32,A33; - A11=0;A12=0;A13=0;A21=0;A22=0;A23=0;A31=0;A32=0;A33=0; + Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33; + a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0; Double_t xbar = 0; Double_t ybar = 0; Double_t zbar = 0; @@ -746,29 +785,29 @@ Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Doubl Double_t value = 0.0; // minimum eigenvalue Double_t x1,x2,x3; // eigenvector component Double_t n1,n2,n3,nr= 0;// unit eigenvector - Double_t Radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm] - Double_t sigma_MS = 0; - TVector3 Vec,VecNor; + Double_t radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm] + Double_t sigmaMS = 0; + TVector3 vVec,vVecNor; // Select RecPoints belonging to the track for(k =0;ktrack==tracknumber) { - name = fPointRecs[k]->fName; - Pt = fPointRecs[k]->fPt; - Layer = fPointRecs[k]->lay; - Int_t ilay = Layer-1; + if(fPointRecs[k]->GetTrack()==tracknumber) { + name = fPointRecs[k]->GetName(); + pt = fPointRecs[k]->GetPt(); + pLayer = fPointRecs[k]->GetLay(); + Int_t ilay = pLayer-1; if(npl[ilay]!=0) continue; if(bitlay[ilay] == 1) { - original[N]->SetXYZ(0.1*fPointRecs[k]->fx,0.1*fPointRecs[k]->fy,0.1*fPointRecs[k]->fz); - errors[N]->SetXYZ(0.1*fPointRecs[k]->fdx,0.1*fPointRecs[k]->fdy,0.1*fPointRecs[k]->fdz); - sigma_MS = (Radiusdm[Layer]-Radiusdm[0])*0.000724/P;// beam pipe contribution - for(Int_t j=1;jSetXYZ(0.1*fPointRecs[k]->GetX(),0.1*fPointRecs[k]->GetY(),0.1*fPointRecs[k]->GetZ()); + errors[nN]->SetXYZ(0.1*fPointRecs[k]->GetdX(),0.1*fPointRecs[k]->GetdY(),0.1*fPointRecs[k]->GetdZ()); + sigmaMS = (radiusdm[pLayer]-radiusdm[0])*0.000724/pP;// beam pipe contribution + for(Int_t j=1;jPerp2() )*( 1+ original[N]->Perp2() )/ - ( errors[N]->Perp2() + sigma_MS*sigma_MS ); - linerrors[N]->SetXYZ(errors[N]->X(),errors[N]->Y(),sqrt(errors[N]->Z()*errors[N]->Z()+sigma_MS*sigma_MS)); - N++; + weight[nN] = ( 1 + original[nN]->Perp2() )*( 1+ original[nN]->Perp2() )/ + ( errors[nN]->Perp2() + sigmaMS*sigmaMS ); + linerrors[nN]->SetXYZ(errors[nN]->X(),errors[nN]->Y(),sqrt(errors[nN]->Z()*errors[nN]->Z()+sigmaMS*sigmaMS)); + nN++; npl[ilay]++; } // end if on layer } //end if track==tracknumber @@ -778,7 +817,7 @@ Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Doubl // if(original[5]->X() == 9999 || original[6]->X() != 9999) { - delete [] Weight; + delete [] weight; return 1; // not enough points } @@ -788,54 +827,54 @@ Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Doubl // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE // - RiemannTransf(N,original,Riemann); + RiemannTransf(nN,original,riemann); - Double_t Sum_Weights = 0.0; // sum of weights factor + Double_t sumWeights = 0.0; // sum of weights factor - for(Int_t j=0;jX(); - ybar+=Weight[j]*Riemann[j]->Y(); - zbar+=Weight[j]*Riemann[j]->Z(); - Sum_Weights+=Weight[j]; + for(Int_t j=0;jX(); + ybar+=weight[j]*riemann[j]->Y(); + zbar+=weight[j]*riemann[j]->Z(); + sumWeights+=weight[j]; } - xbar /= Sum_Weights; - ybar /= Sum_Weights; - zbar /= Sum_Weights; + xbar /= sumWeights; + ybar /= sumWeights; + zbar /= sumWeights; - for(Int_t j=0;jX() - xbar)*(Riemann[j]->X() - xbar); - A12 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->Y() - ybar); - A22 += Weight[j]*(Riemann[j]->Y() - ybar)*(Riemann[j]->Y() - ybar); - A23 += Weight[j]*(Riemann[j]->Y() - ybar)*(Riemann[j]->Z() - zbar); - A13 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->Z() - zbar); - A33 += Weight[j]*(Riemann[j]->Z() - zbar)*(Riemann[j]->Z() - zbar); + for(Int_t j=0;jX() - xbar)*(riemann[j]->X() - xbar); + a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar); + a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar); + a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar); + a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar); + a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar); } - A11 /= N; - A12 /= N; - A22 /= N; - A23 /= N; - A13 /= N; - A33 /= N; - A21 = A12; - A32 = A23; - A31 = A13; + a11 /= nN; + a12 /= nN; + a22 /= nN; + a23 /= nN; + a13 /= nN; + a33 /= nN; + a21 = a12; + a32 = a23; + a31 = a13; // ************** Determinant parameters ******************** // n.b. simplifications done keeping in mind symmetry of A // a = 1; - b = (-A11-A33-A22); - c = (A11*(A22+A33)+A33*A22-A12*A21-A13*A31-A23*A32); - d = (A31*A22*A13+(A12*A21-A11*A22)*A33-2.0*A23*A13*A12+A11*A23*A32); + b = (-a11-a33-a22); + c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32); + d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32); // ************** Find the 3 eigenvalues ************************* - Int_t Check_Cubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]); + Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]); - if(Check_Cubic !=1 ){ + if(checkCubic !=1 ){ printf("Track %d Has no real solution continuing ...\n",tracknumber); - delete [] Weight; + delete [] weight; return 2; } @@ -846,13 +885,13 @@ Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Doubl // ************ Eigenvector relative to value ************** x3 = 1; - x2 = (A33*A21-A23*A31-value*A21)/(A22*A31-A32*A21-value*A31); - x1 = (value-A33-A32*x2)/A31; - Vec.SetXYZ(x1,x2,x3); - VecNor = Vec.Unit(); - n1 = VecNor.X(); - n2 = VecNor.Y(); - n3 = VecNor.Z(); + x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31); + x1 = (value-a33-a32*x2)/a31; + vVec.SetXYZ(x1,x2,x3); + vVecNor = vVec.Unit(); + n1 = vVecNor.X(); + n2 = vVecNor.Y(); + n3 = vVecNor.Z(); nr = -n1*xbar-n2*ybar-n3*zbar; u0 = -0.5*n1/(nr+n3); @@ -885,212 +924,379 @@ Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Doubl rho *= 10.0; u0 *= 10.0; v0 *= 10.0; - ierrl=FitLinear(N,original,linerrors,omega,u0,v0,fphi,zData,zError,CorrLin); + ierrl=FitLinear(nN,original,linerrors,omega,u0,v0,fphi,zData,zError,corrLin); chisql=zError.Z(); // fprintf(pout,"%f \n",chisql); z0=zData.X(); vpar=zData.Y(); - fCorrLin = CorrLin; + fCorrLin = corrLin; ierr = (ierrl > ierr ? ierrl : ierr); // fclose(pout); - delete [] Weight; + delete [] weight; return ierr; } +Int_t AliITSRiemannFit::FitHelix(Int_t NPoints, TVector3** fPointRecs,TVector3** fPointRecErrors,Float_t& f1, Float_t& f2, Float_t& f3) { -//----------------------------------------------------------------------------- -void AliITSRiemannFit::Streamer(TBuffer &lRb){ -//////////////////////////////////////////////////////////////////////// -// The default Streamer function "written by ROOT" doesn't write out -// the arrays referenced by pointers. Therefore, a specific Streamer function -// has to be written. This function should not be modified but instead added -// on to so that older versions can still be read. -//////////////////////////////////////////////////////////////////////// - // Stream an object of class AliITSRiemannFit. - Int_t i,j,n; - Int_t ii,jj; - n=20; - - if (lRb.IsReading()) { - Version_t lRv = lRb.ReadVersion(); if (lRv) { } - TObject::Streamer(lRb); - lRb >> fSizeEvent; - lRb >> fPrimaryTracks; - lRb >> fPoints; - for(i=0;i<6;i++) lRb >> fPLay[i]; - if(fPointRecs!=0){ - for(i=0;i> fPointRecs[i][j].lay; - lRb >> fPointRecs[i][j].lad; - lRb >> fPointRecs[i][j].det; - lRb >> fPointRecs[i][j].track; - lRb >> fPointRecs[i][j].fx; - lRb >> fPointRecs[i][j].fy; - lRb >> fPointRecs[i][j].fz; - lRb >> fPointRecs[i][j].fr; - lRb >> fPointRecs[i][j].fdE; - lRb >> fPointRecs[i][j].fdx; - lRb >> fPointRecs[i][j].fdy; - lRb >> fPointRecs[i][j].fdz; - lRb >> fPointRecs[i][j].fPt; - lRb >> (Char_t*)fPointRecs[i][j].fName; - for (ii=0;ii<4;ii++) - lRb << fPointRecs[i][j].fOrigin[ii]; - for (jj=0;jj<4;jj++) - lRb << fPointRecs[i][j].fMomentum[jj]; - lRb >> fPointRecs[i][j].fCode; - lRb >> fPointRecs[i][j].phi; - lRb >> fPointRecs[i][j].eta; - lRb >> fPointRecs[i][j].vertexPhi; - } //end for j - } //end for i -// if(fspdi!=0){ -// for(i=0;i> fspdi[i][j].lay; -// lRb >> fspdi[i][j].lad; -// lRb >> fspdi[i][j].det; -// lRb >> fspdi[i][j].track; -// lRb >> fspdi[i][j].fx; -// lRb >> fspdi[i][j].fy; -// lRb >> fspdi[i][j].fz; -// lRb >> fspdi[i][j].fr; -// lRb >> fspdi[i][j].fdE; -// lRb >> fspdi[i][j].fdx; -// lRb >> fspdi[i][j].fdy; -// lRb >> fspdi[i][j].fdz; -// lRb >> fspdi[i][j].fPt; -// for (ii=0;ii<4;ii++) -// lRb << fspdi[i][j].fOrigin[ii]; -// for (jj=0;jj<4;jj++) -// lRb << fspdi[i][j].fMomentum[jj]; -// lRb >> fspdi[i][j].fCode; -// lRb >> (Char_t*)fspdi[i][j].fName; -// lRb >> fspdi[i][j].phi; -// lRb >> fspdi[i][j].eta; -// lRb >> fspdi[i][j].vertexPhi; -// } //end for j -// } //end for i -// if(fspdo!=0){ -// for(i=0;i> fspdo[i][j].lay; -// lRb >> fspdo[i][j].lad; -// lRb >> fspdo[i][j].det; -// lRb >> fspdo[i][j].track; -// lRb >> fspdo[i][j].fx; -// lRb >> fspdo[i][j].fy; -// lRb >> fspdo[i][j].fz; -// lRb >> fspdo[i][j].fr; -// lRb >> fspdo[i][j].fdE; -// lRb >> fspdo[i][j].fdx; -// lRb >> fspdo[i][j].fdy; -// lRb >> fspdo[i][j].fdz; -// lRb >> fspdo[i][j].fPt; -// for (ii=0;ii<4;ii++) -// lRb << fspdo[i][j].fOrigin[ii]; -// for (jj=0;jj<4;jj++) -// lRb << fspdo[i][j].fMomentum[jj]; -// lRb >> fspdo[i][j].fCode; -// lRb >> (Char_t*)fspdo[i][j].fName; -// lRb >> fspdo[i][j].phi; -// lRb >> fspdo[i][j].eta; -// lRb >> fspdo[i][j].vertexPhi; -// } //end for j -// } //end for i - } else { - lRb.WriteVersion(AliITSRiemannFit::IsA()); - TObject::Streamer(lRb); - lRb << fSizeEvent; - lRb << fPrimaryTracks; - lRb << fPoints; - for(i=0;i<6;i++) lRb >> fPLay[i]; - for(i=0;iSetXYZ(9999,9999,9999); + } + + // + // Riemann fit parameters + // + Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33; + a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0; + Double_t xbar = 0; + Double_t ybar = 0; + Double_t zbar = 0; + // + Double_t a,b,c,d; // cubic parameters + Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions + Double_t value = 0.0; // minimum eigenvalue + Double_t x1,x2,x3; // eigenvector component + Double_t n1,n2,n3,nr= 0; // unit eigenvector + TVector3 vVec,vVecNor; + + for (Int_t ip=0; ipSetXYZ(0.1*fPointRecs[ip]->X(),0.1*fPointRecs[ip]->Y(),0.1*fPointRecs[ip]->Z()); + +errors[ip]->SetXYZ(0.1*fPointRecErrors[ip]->X(),0.1*fPointRecErrors[ip]->Y(),0.1*fPointRecErrors[ip]->Z()); + weight[ip] = (1+original[ip]->Perp2())*(1+original[ip]->Perp2())/(errors[ip]->Perp2()); + linerrors[ip]->SetXYZ(errors[ip]->X(),errors[ip]->Y(),errors[ip]->Z()); + } + + + // + // + // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE + // + + RiemannTransf(NPoints,original,riemann); + + Double_t sumWeights = 0.0; // sum of weights factor + + for(Int_t j=0;jX(); + ybar+=weight[j]*riemann[j]->Y(); + zbar+=weight[j]*riemann[j]->Z(); + sumWeights+=weight[j]; + } + + xbar /= sumWeights; + ybar /= sumWeights; + zbar /= sumWeights; + + for(Int_t j=0;jX() - xbar)*(riemann[j]->X() - xbar); + a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar); + a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar); + a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar); + a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar); + a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar); + } + // + // this doesn't seem to work... + // + // a11 /= sumWeights; + // a12 /= sumWeights; + // a22 /= sumWeights; + // a23 /= sumWeights; + // a13 /= sumWeights; + // a33 /= sumWeights; + + a11 /= NPoints; + a12 /= NPoints; + a22 /= NPoints; + a23 /= NPoints; + a13 /= NPoints; + a33 /= NPoints; + a21 = a12; + a32 = a23; + a31 = a13; + +// ************** Determinant parameters ******************** +// n.b. simplifications done keeping in mind symmetry of A +// + a = 1; + b = (-a11-a33-a22); + c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32); + d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32); + +// ************** Find the 3 eigenvalues ************************* + Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]); + + if(checkCubic !=1 ){ + printf("No real solution. Check data.\n"); + delete [] weight; + return 999; + } + +// **************** Find the lowest eigenvalue ***************** + if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0]; + if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1]; + if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2]; + + // ************ Eigenvector relative to value ************** + x3 = 1; + x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31); + x1 = (value-a33-a32*x2)/a31; + vVec.SetXYZ(x1,x2,x3); + vVecNor = vVec.Unit(); + n1 = vVecNor.X(); + n2 = vVecNor.Y(); + n3 = vVecNor.Z(); + nr = -n1*xbar-n2*ybar-n3*zbar; + + u0 = -0.5*n1/(nr+n3); + v0 = -0.5*n2/(nr+n3); + rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3))); + + + fFit = 0.0; + for (Int_t i=0; iX()-u0)*(original[i]->X()-u0)+(original[i]->Y()-v0)*(original[i]->Y()-v0))-rho); + } + fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns + fphi = TMath::ATan2(v0,u0); + +//************************************************************************** +// LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s +// strictly linear (no approximation) +//************************************************************************** + + //////////////////////////////////////////////////////////////////////////////////////////////////////////// + // // + // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S // + // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... // + // // + //////////////////////////////////////////////////////////////////////////////////////////////////////////// + + rho *= 10.0; + u0 *= 10.0; + v0 *= 10.0; + + ierrl=LinearFit(NPoints,original,linerrors,kOmega,u0,v0,fphi,zData,zError,corrLin); + if(ierrl==33) return 0; + chisql=zError.Z(); +// fprintf(pout,"%f \n",chisql); + z0=zData.X(); + vpar=zData.Y(); + fCorrLin = corrLin; + ierr = (ierrl > ierr ? ierrl : ierr); +// fclose(pout); + delete [] weight; + + f1=fphi; + f2=vpar/(kOmega*TMath::Abs(rho)); + f3=1/rho; + delete[] ori; + delete[] rie; + delete[] err; + delete[] linerr; + delete[] original; + delete[] riemann; + delete[] errors; + delete[] linerrors; + + return 1; + + +} + +//____________________________________________________________ + +Int_t AliITSRiemannFit::LinearFit(Int_t npoints, TVector3 **input, + TVector3 **errors, Double_t omega, + Double_t &thu0, Double_t &thv0, Double_t &phi,TVector2 &zData, TVector3 &zError, + Double_t &corrLin){ + /////////////////////////////////////////////////////////////////////// + // Fit the points in the (z,s) plane - helix 3rd equation + // + /////////////////////////////////////////////////////////////////////// + //By R.Turrisi + + Int_t direction=0; + //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints]; + //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints]; + Double_t * z = new Double_t[npoints]; + Double_t * x = new Double_t[npoints]; + Double_t * y = new Double_t[npoints]; + Double_t * s = new Double_t[npoints]; + Double_t * ez = new Double_t[npoints]; + Double_t * ex = new Double_t[npoints]; + Double_t * ey = new Double_t[npoints]; + Double_t * es = new Double_t[npoints]; + Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare; + + + // Double_t chi=TMath::Pi()/2.0+phi; + Double_t chi=-TMath::Pi()-phi; + Double_t angold=0.0, tpang=0.0; + for(Int_t k = 0; kX(); ex[k] = 10.0*errors[k]->X(); + y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y(); + z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z(); + if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles... + chisquare=9999.99; + cerr<<"limit for x-x_0 "<0) direction+=(z[k]>z[k-1] ? 1 : -1); + s[k] = (ang1+chi)/omega; + es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]); + } + if ( TMath::Abs(direction) != (npoints-1) ) {return 11;} + + // if(s[0]>-636 && s[0]<-625) return 33; + + TGraph* fitHist = new TGraph(npoints,s,z); + TF1* f1 = new TF1("f1",Fitfunction,-100,100,2); + + f1->SetParameter(0,1); + f1->SetParameter(1,1); + f1->SetLineColor(2); + fitHist->Fit(f1,"qQ"); + + z0 = f1->GetParameter(0); + vpar = f1->GetParameter(1); + ez0 = f1->GetParError(0); + evpar= f1->GetParError(1); + chisquare=f1->GetChisquare(); + zData.Set(z0,vpar); + zError.SetXYZ(ez0,evpar,chisquare); + + Double_t sigmas=0.; + Double_t sigmaz=0.; + Double_t avs=0.; + Double_t avz=0.; + Double_t avsz=0.; + + for(Int_t j = 0; j < npoints; j++) { + avs += s[j]; + avz += z[j]; + avsz += s[j]*z[j]; + } + avs /= (Double_t)npoints; + avz /= (Double_t)npoints; + avsz /= (Double_t)npoints; + + for(Int_t l = 0; l < npoints; l++) { + sigmas += (s[l]-avs)*(s[l]-avs); + sigmaz += (z[l]-avz)*(z[l]-avz); + } + sigmas /=(Double_t)npoints; + sigmaz /=(Double_t)npoints; + + sigmas = sqrt(sigmas); + sigmaz = sqrt(sigmaz); + + corrLin = (avsz-avs*avz)/(sigmas*sigmaz); + + + + delete [] z; + delete [] x; + delete [] y; + delete [] s; + delete [] ez; + delete [] ex; + delete [] ey; + delete [] es; + delete f1; delete fitHist; + return 0; } + + +//_______________________________________________________ + +Double_t AliITSRiemannFit::Fitfunction(Double_t *x, Double_t* par){ + // function used for fit + return par[0]+(*x)*par[1]; + +} +