1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
16 * ///////////////////////////////////////////////////////////////////// *
18 * This class performs a fast fit of helices going through the <=6 *
19 * points of the ITS, with the goal of studying tracking and *
20 * vertexing performances. *
21 * Generated kinematics is used to take into account different weights *
22 * associated to points in different layers (with different multiple *
23 * scattering-originated errors). *
25 * Based on the work by A. Strandlie, R. Fruhwirth *
27 * First implementation by N. Bustreo, R. Turrisi - July 2000 *
29 * Further modifications by A. Dainese, R. Turrisi *
31 * Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it *
33 * **************************************************************************/
36 // Modified November, 7th 2001 by Rosario Turrisi
37 // (rosario.turrisi@pd.infn.it)
39 // FitHelix returns different values. 0=ok, >0 =problem
40 // void FitLinear -> Int_t FitLinear to give feedback of errors to FitHelix
43 // Modified July, 30th 2001 by Rosario Turrisi
44 // (rosario.turrisi@pd.infn.it)
46 // Fit for z now in (z,s) plane.
47 // Returns parameters in order to write the helix equation
48 // and find the right phase/initial point.
50 // "PROPER WEIGHTS": (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2)
53 #include "AliITSRiemannFit.h"
55 #include "TClonesArray.h"
61 #include "TGraphErrors.h"
66 #include "TParticle.h"
68 #include "AliITSRecPoint.h"
69 #include "AliITSgeom.h"
70 #include "AliITSmodule.h"
72 ClassImp(AliITSRiemannFit)
75 AliITSRiemannFit::AliITSRiemannFit() {
76 ///////////////////////////////////////////////////////////
77 // Default constructor.
78 // Set everything to zero.
79 ////////////////////////////////////////////////////////////
89 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
92 //----------------------------------------------------------------------
94 AliITSRiemannFit::~AliITSRiemannFit() {
95 ///////////////////////////////////////////////////////////
96 // Default destructor.
97 // if arrays exist delete them. Then set everything to zero.
98 ////////////////////////////////////////////////////////////
100 for(Int_t i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
102 } // end if fPointRecs!=0
111 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
114 //----------------------------------------------------------------------
116 AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks) {
117 ///////////////////////////////////////////////////////////
119 // Set fSizeEvent to size and fPrimaryTracks to ntracks.
121 ////////////////////////////////////////////////////////////
125 fPrimaryTracks = ntracks;
130 Point_tl *first = new Point_tl[fSizeEvent];
131 Point_tl **PointRecs = new Point_tl*[fSizeEvent];
132 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
133 for(Int_t j=0;j<fSizeEvent;j++) // create an array of struct
134 PointRecs[j] = &(first[j]);
136 // ---------------------------------------------------------------------
138 void FillPoints(Point_tl **Points,Int_t &index,Float_t *xpoint,
140 TLorentzVector PE,TLorentzVector OT,Int_t *id,
141 Int_t track,const Char_t *name,Int_t code,
143 ///////////////////////////////////////////////////////////////////////
144 // Fill the structure Point_tl with the proper data
146 //////////////////////////////////////////////////////////////////////
147 Float_t PI2 = 2.0*TMath::Pi();
155 phi = TMath::ATan2(y,x);
156 if(phi<0.0) phi += PI2;
157 Points[i]->phi = phi;
158 Points[i]->eta = -0.5*tan(0.5*TMath::ATan2(r,z));
162 Points[i]->fdx = error[0];
163 Points[i]->fdy = error[1];
164 Points[i]->fdz = error[2];
166 Points[i]->track = track;
167 Points[i]->lay = id[0],
168 Points[i]->lad = id[1];
169 Points[i]->det = id[2];
170 Points[i]->fMomentum = PE;
171 Points[i]->fOrigin = OT;
172 Points[i]->fPt = sqrt(PE.X()*PE.X()+PE.Y()*PE.Y());
173 Points[i]->fCode = code;
174 Points[i]->fName = name;
175 Points[i]->vertexPhi = phiorigin;
180 // -----------------------------------------------------------------------
182 void AliITSRiemannFit::InitPoints(Int_t evnt,Int_t ntracks,AliITS *ITS,
183 TTree *TR,Int_t nparticles){
184 //////////////////////////////////////////////////////////////////////
185 // Fill the class member fPointRecs with the reconstructed points
186 // Set All other members to the real values
188 /////////////////////////////////////////////////////////////////////
189 printf("\n ************* Starting Init Points *************\n");
191 AliITSgeom *gm = (AliITSgeom*)ITS->GetITSgeom();
192 //get pointer to modules array
193 TObjArray *ITSmodules = ITS->GetModules();
194 Int_t nmodules=ITSmodules->GetEntriesFast();
195 printf("nmodules = %d \n",nmodules);
196 // Get the points from points file
197 AliITSmodule *itsModule;
200 AliITSRecPoint *recp;
201 nent=TR->GetEntries();
202 TClonesArray *ITSrec = ITS->RecPoints();
205 for (mod=0; mod<nmodules; mod++) {
206 itsModule=(AliITSmodule*)ITSmodules->At(mod);
207 ITS->ResetRecPoints();
209 Int_t nrecp = ITSrec->GetEntries();
215 fPrimaryTracks = ntracks;
216 fParticles = nparticles;
217 Point_tl *global = new Point_tl[iMAX];
218 fPointRecs = new Point_tl*[iMAX];
221 // Point_tl *first = new Point_tl[iMAX];
222 // Point_tl *second = new Point_tl[iMAX];
223 // fspdi = new Point_tl*[iMAX];
224 // fspdo = new Point_tl*[iMAX];
225 for(Int_t j=0;j<iMAX;j++) {
226 fPointRecs[j] = &(global[j]);
229 // fspdi[j] = &(first[j]);
230 // fspdo[j] = &(second[j]);
233 Int_t ieta=0,ieta2=0;
234 Int_t i,id[4],idold[4];
235 Int_t track=0;// // track of hit
236 Float_t xpoint[3],error_plus[3],error_minus[3],global_error[3]; // position and error of the point
237 TLorentzVector OT,PE;
238 Float_t locals[3],locals_error[3],locals_plus[3],locals_minus[3]; // local position and local errors
242 Int_t layer,ladder,detector;
243 Float_t xcluster,zcluster;
244 Int_t num=0,nspdi=0,nspdo=0,nsddi=0,nsddo=0,nssdi=0,nssdo=0;
246 for (mod=0; mod<nmodules; mod++) {
247 itsModule=(AliITSmodule*)ITSmodules->At(mod);
248 ITS->ResetRecPoints();
250 Int_t nrecp = ITSrec->GetEntries();
251 if (!nrecp) continue;
252 itsModule->GetID(layer,ladder,detector);
254 for (irec=0;irec<nrecp;irec++) {
255 recp = (AliITSRecPoint*)ITSrec->UncheckedAt(irec);
256 track=recp->fTracks[0];
257 if(track <0 ) continue;
258 xcluster=recp->GetX(); // x on cluster
259 zcluster=recp->GetZ(); // z on cluster
260 part = (TParticle*) gAlice->Particle(track);
261 part->ProductionVertex(OT); // set the vertex
262 part->Momentum(PE); // set the vertex momentum
263 name = part->GetName();
264 code = part->GetPdgCode();
270 locals[0]=xcluster; // x on cluster
271 locals[1]=0.0; // y on cluster
272 locals[2]=zcluster; // z on cluster
273 locals_error[0]=sqrt(recp->GetSigmaX2());
275 locals_error[2]=sqrt(recp->GetSigmaZ2());
276 locals_plus[0]=xcluster+sqrt(recp->GetSigmaX2()); // x on cluster
277 if(layer==1||layer==2) locals_plus[1]=0.0150/2; // y on cluster
278 else if(layer==3||layer==4) locals_plus[1]=0.0280/2; // y on cluster
279 else if(layer==5||layer==6) locals_plus[1]=0.0300/2; // y on cluster
280 locals_plus[2]=zcluster+sqrt(recp->GetSigmaZ2()); // z on cluster
281 locals_minus[0]=xcluster-sqrt(recp->GetSigmaX2()); // x on cluster
282 locals_minus[1]=0.0; // y on cluster
283 locals_minus[2]=zcluster-sqrt(recp->GetSigmaZ2()); // z on cluster
285 gm->LtoG(layer,ladder,detector,locals,xpoint);
286 gm->LtoG(layer,ladder,detector,locals_plus,error_plus);
287 gm->LtoG(layer,ladder,detector,locals_minus,error_minus);
288 global_error[0]=0.5*TMath::Abs(error_plus[0]-error_minus[0]);
289 global_error[1]=0.5*TMath::Abs(error_plus[1]-error_minus[1]);
290 global_error[2]=0.5*TMath::Abs(error_plus[2]-error_minus[2]);
292 if(TMath::Abs(part->Eta())<=1.0) ieta++;
293 if(TMath::Abs(part->Eta())<=0.5) ieta2++;
295 if(!(id[0]==idold[0]&&id[1]==idold[1]&&
296 id[2]==idold[2]&&id[3]==idold[3])) {
297 FillPoints(fPointRecs,num,xpoint,global_error,PE,OT,id,track,name,code,Phi);
321 // FillPoints(fspdi,nspdi,xpoint,global_error,PE,OT,id,track,name,code,Phi);
325 // FillPoints(fspdo,nspdo,xpoint,global_error,PE,OT,id,track,name,code,Phi);
339 for(i=0;i<4;i++) idold[i] = id[i];
340 for(i=0;i<3;i++) xpoint[i] = 0.0;
341 } // end if id != idold
353 printf("%d primary tracks in eta=+-1\n",ieta);
354 printf("%d primary tracks#2 in eta=+-0.5\n",ieta2);
355 printf("\nInitPoints :\n\nPoints on Layer1 : %d on Layer2 : %d\n",nspdi,nspdo);
356 printf("Points on Layer3 : %d on Layer4 : %d\n",nsddi,nsddo);
357 printf("Points on Layer5 : %d on Layer6 : %d\n",nssdi,nssdo);
358 printf("Points on all Layers: %d\n",num);
359 printf("\n ************* Init Points Finished *************\n");
362 // ------------------------------------------------------------------------
363 ///////////////////////////////////////////////////////////
364 // Functions for sorting the fPointRecs array
365 ///////////////////////////////////////////////////////////
366 Bool_t SortZ(const Point_tl *s1,const Point_tl *s2){
367 // Z sorting function for qsort.
371 if(a<0.0) return kTRUE;
372 if(a>0.0) return kFALSE;
375 Bool_t SortTrack(const Point_tl *s1,const Point_tl *s2){
376 // track sorting function for qsort.
379 a = s1->track - s2->track;
380 if(a<0.0) return kTRUE;
381 if(a>0.0) return kFALSE;
384 void hpsortTrack(Point_tl **ra,Int_t n){
390 l = ((n-1) >> 1) +1; // divide 2 + 1
394 rra = ra[--l]; // decrement first
398 if(--ir == 0){ // decrement first
406 if( j<ir && SortTrack(ra[j],ra[j+1]) ) j++;
407 if( SortTrack(rra,ra[j]) ){
418 void hpsortZ(Point_tl **ra,Int_t n){
424 l = ((n-1) >> 1) +1; // devide 2 + 1
428 rra = ra[--l]; // decrament first
432 if(--ir == 0){ // decrament first
440 if( j<ir && SortZ(ra[j],ra[j+1]) ) j++;
441 if( SortZ(rra,ra[j]) ){
452 //-----------------------------------------------------------------------
453 ////////////////////////////////////////////////////////////////////
455 ///////////////////////////////////////////////////////////////////
456 Int_t Partition(Int_t array[],Int_t left,Int_t right){
457 Int_t val = array[left];
459 Int_t rm = right + 1;
470 Int_t tempr = array[rm];
481 ///////////////////////////////////////////////////////////////////////
483 void AliITSRiemannFit::WritePoints(void) {
484 /////////////////////////////////////////////////////////////////////
485 // write the data in a file (temporary ascii)
486 /////////////////////////////////////////////////////////////////////
487 FILE *ascii= fopen("AsciiPoints.dat","w");
488 for(Int_t i=0;i<fPoints;i++) {
489 fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->lay,
490 fPointRecs[i]->track,fPointRecs[i]->fx,fPointRecs[i]->fy,
496 //-----------------------------------------------------------------------
498 void AliITSRiemannFit::ReadPoints(void) {
499 //////////////////////////////////////////////////////////////////////
500 // read the filled array
501 /////////////////////////////////////////////////////////////////////
502 hpsortTrack(fPointRecs,fPoints);
503 for(Int_t i=0;i<fPoints;i++)
504 printf("%d\t%d\t%d\t%f\t%f\t%f\t(%.0f,%.0f,%.0f)\t%.3f\t%s\n",
505 i,fPointRecs[i]->lay,fPointRecs[i]->track,fPointRecs[i]->fx,
506 fPointRecs[i]->fy,fPointRecs[i]->fz,fPointRecs[i]->fOrigin.X(),
507 fPointRecs[i]->fOrigin.Y(),fPointRecs[i]->fOrigin.Z(),
508 fPointRecs[i]->fPt,fPointRecs[i]->fName);
511 //-----------------------------------------------------------------------
513 Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c,
514 Double_t &x1,Double_t &x2,Double_t &x3){
515 //////////////////////////////////////////////
516 /// Solve cubic equation:
517 /// x^3 + a*x^2 +b*x + c
519 /// returns x1 , x2 , x3
520 ////////////////////////////////////////
522 Double_t Q = ((a*a - 3*b)/9);
523 Double_t R = ((2*a*a*a - 9*a*b +27*c)/54);
525 Double_t F = -2*sqrt(Q);
527 Double_t PI2 = TMath::Pi()*2;
530 cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<<endl;
537 theta = TMath::ACos( R/sqrt(Q*Q*Q));
539 x1 = (F*TMath::Cos(theta/3))-g;
540 x2 = (F*TMath::Cos((theta+PI2)/3))-g;
541 x3 = (F*TMath::Cos((theta-PI2)/3))-g;
545 //-----------------------------------------------------------------
547 void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) {
548 ///////////////////////////////////////////////////////////////////////
549 // This function apllies the transformation in the Riemann sphere
551 ///////////////////////////////////////////////////////////////////////
552 Float_t *R = new Float_t[npoints];
553 Float_t *Theta = new Float_t[npoints];
554 Float_t PI2 = 2*TMath::Pi();
557 for(Int_t i=0;i<npoints;i++) {
558 R[i] = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y());
559 Theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X());
560 if(Theta[i]<0) Theta[i]+=PI2;
561 x = R[i]*cos(Theta[i])/(1+R[i]*R[i]);
562 y = R[i]*sin(Theta[i])/(1+R[i]*R[i]);
563 z = R[i]*R[i]/(1+R[i]*R[i]);
564 To[i]->SetXYZ(x,y,z);
572 //---------------------------------------------------------------------
574 Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t omega,
575 Double_t &thu0, Double_t &thv0, Double_t &phi, TVector2 &zData, TVector3 &zError,
577 ///////////////////////////////////////////////////////////////////////
578 // Fit the points in the (z,s) plane - helix 3rd equation
580 ///////////////////////////////////////////////////////////////////////
582 //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
583 //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
584 Double_t * z = new Double_t[npoints];
585 Double_t * x = new Double_t[npoints];
586 Double_t * y = new Double_t[npoints];
587 Double_t * s = new Double_t[npoints];
588 Double_t * ez = new Double_t[npoints];
589 Double_t * ex = new Double_t[npoints];
590 Double_t * ey = new Double_t[npoints];
591 Double_t * es = new Double_t[npoints];
592 Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
594 // Double_t chi=TMath::Pi()/2.0+phi;
595 Double_t chi=-TMath::Pi()-phi;
596 Double_t angold=0.0, tpang=0.0;
597 for(Int_t k = 0; k<npoints; k++) {
598 x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
599 y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
600 z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
601 if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
603 cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
614 Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
615 if( (x[k]-thu0)<0 ) {
617 tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
622 if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
623 s[k] = (ang1+chi)/omega;
624 es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
626 if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
628 TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez);
629 fitHist->Fit("pol1","Q");
630 z0 = fitHist->GetFunction("pol1")->GetParameter(0);
631 vpar = fitHist->GetFunction("pol1")->GetParameter(1);
632 ez0 = fitHist->GetFunction("pol1")->GetParError(0);
633 evpar = fitHist->GetFunction("pol1")->GetParError(1);
634 chisquare = fitHist->GetFunction("pol1")->GetChisquare();
636 zError.SetXYZ(ez0,evpar,chisquare);
644 for(Int_t j = 0; j < npoints; j++) {
649 Avs /= (Double_t)npoints;
650 Avz /= (Double_t)npoints;
651 Avsz /= (Double_t)npoints;
653 for(Int_t l = 0; l < npoints; l++) {
654 Sigmas += (s[l]-Avs)*(s[l]-Avs);
655 Sigmaz += (z[l]-Avz)*(z[l]-Avz);
657 Sigmas /=(Double_t)npoints;
658 Sigmaz /=(Double_t)npoints;
660 Sigmas = sqrt(Sigmas);
661 Sigmaz = sqrt(Sigmaz);
663 CorrLin = (Avsz-Avs*Avz)/(Sigmas*Sigmaz);
677 //-------------------------------------------------------------------
678 Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Int_t charge,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
679 Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,Double_t& omega, Double_t& z0,
680 Double_t& vpar,Double_t& chisql, Double_t& fCorrLin,Double_t& fFit,
681 Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
682 ///////////////////////////////////////////////////////////////////////
683 // This function finds the helix paramenters
684 // d0 = impact parameter
685 // rho = radius of circle
688 // starting from the momentum and the outcome of
689 // the fit on the Riemann sphere (i.e. u0,v0,rho)
691 // MIND !!!! Here we assume both angular velocities be 1.0 (yes, one-dot-zero !)
694 ///////////////////////////////////////////////////////////////////////
696 // All this stuff relies on this hypothesis !!!
698 // FILE *pout=fopen("chisql.dat","a");
699 Int_t ierr = 0, ierrl=0;
702 Int_t bitlay[6]={1,1,1,1,1,1};
703 bitlay[0]*=first; bitlay[1]*=second; bitlay[2]*=third; bitlay[3]*=fourth; bitlay[4]*=fifth; bitlay[5]*=sixth;
704 fd0 = -9999; // No phisycs value
705 u0 = -9999.9999; // parameters of helix - strange value...
706 v0 = -9999.9999; // parameters of helix - strange value...
707 rho = -9999.9999; // parameters of helix -unphysical strange value...
709 const Char_t* name = 0;
713 Int_t npl[6]={0,0,0,0,0,0};
714 Double_t P = sqrt(Px*Px+Py*Py+Pz*Pz);
715 Double_t Pt = sqrt(Px*Px+Py*Py);
719 TVector3 *ori = new TVector3[iMAX];
720 TVector3 **original = new TVector3*[iMAX];
721 TVector3 *rie = new TVector3[iMAX];
722 TVector3 **Riemann = new TVector3*[iMAX];
723 TVector3 *err = new TVector3[iMAX];
724 TVector3 **errors = new TVector3*[iMAX];
725 TVector3 *linerr = new TVector3[iMAX];
726 TVector3 **linerrors = new TVector3*[iMAX];
727 //PH Double_t Weight[iMAX];
728 Double_t * Weight = new Double_t[iMAX];
731 original[i] = &(ori[i]);
732 Riemann[i] = &(rie[i]);
733 errors[i] = &(err[i]);
734 linerrors[i] = &(linerr[i]);
736 for(k =0;k<iMAX;k++) original[k]->SetXYZ(9999,9999,9999);
737 Double_t A11,A12,A13,A21,A22,A23,A31,A32,A33;
738 A11=0;A12=0;A13=0;A21=0;A22=0;A23=0;A31=0;A32=0;A33=0;
742 Double_t a,b,c,d; // cubic parameters
743 Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
744 Double_t value = 0.0; // minimum eigenvalue
745 Double_t x1,x2,x3; // eigenvector component
746 Double_t n1,n2,n3,nr= 0;// unit eigenvector
747 Double_t Radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm]
748 Double_t sigma_MS = 0;
751 // Select RecPoints belonging to the track
752 for(k =0;k<fPoints;k++){
753 if(fPointRecs[k]->track==tracknumber) {
754 name = fPointRecs[k]->fName;
755 Pt = fPointRecs[k]->fPt;
756 Layer = fPointRecs[k]->lay;
757 Int_t ilay = Layer-1;
758 if(npl[ilay]!=0) continue;
759 if(bitlay[ilay] == 1) {
760 original[N]->SetXYZ(0.1*fPointRecs[k]->fx,0.1*fPointRecs[k]->fy,0.1*fPointRecs[k]->fz);
761 errors[N]->SetXYZ(0.1*fPointRecs[k]->fdx,0.1*fPointRecs[k]->fdy,0.1*fPointRecs[k]->fdz);
762 sigma_MS = (Radiusdm[Layer]-Radiusdm[0])*0.000724/P;// beam pipe contribution
763 for(Int_t j=1;j<Layer;j++) {
764 sigma_MS += (Radiusdm[Layer]-Radiusdm[j])*0.00136/P;
766 Weight[N] = ( 1 + original[N]->Perp2() )*( 1+ original[N]->Perp2() )/
767 ( errors[N]->Perp2() + sigma_MS*sigma_MS );
768 linerrors[N]->SetXYZ(errors[N]->X(),errors[N]->Y(),sqrt(errors[N]->Z()*errors[N]->Z()+sigma_MS*sigma_MS));
772 } //end if track==tracknumber
775 // 6 points, no more, no less
777 if(original[5]->X() == 9999 || original[6]->X() != 9999)
780 return 1; // not enough points
786 // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
789 RiemannTransf(N,original,Riemann);
791 Double_t Sum_Weights = 0.0; // sum of weights factor
793 for(Int_t j=0;j<N;j++){ // mean values for x[i],y[i],z[i]
794 xbar+=Weight[j]*Riemann[j]->X();
795 ybar+=Weight[j]*Riemann[j]->Y();
796 zbar+=Weight[j]*Riemann[j]->Z();
797 Sum_Weights+=Weight[j];
804 for(Int_t j=0;j<N;j++) { // Calculate the matrix elements
805 A11 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->X() - xbar);
806 A12 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->Y() - ybar);
807 A22 += Weight[j]*(Riemann[j]->Y() - ybar)*(Riemann[j]->Y() - ybar);
808 A23 += Weight[j]*(Riemann[j]->Y() - ybar)*(Riemann[j]->Z() - zbar);
809 A13 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->Z() - zbar);
810 A33 += Weight[j]*(Riemann[j]->Z() - zbar)*(Riemann[j]->Z() - zbar);
823 // ************** Determinant parameters ********************
824 // n.b. simplifications done keeping in mind symmetry of A
828 c = (A11*(A22+A33)+A33*A22-A12*A21-A13*A31-A23*A32);
829 d = (A31*A22*A13+(A12*A21-A11*A22)*A33-2.0*A23*A13*A12+A11*A23*A32);
831 // ************** Find the 3 eigenvalues *************************
832 Int_t Check_Cubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
834 if(Check_Cubic !=1 ){
835 printf("Track %d Has no real solution continuing ...\n",tracknumber);
840 // **************** Find the lowest eigenvalue *****************
841 if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
842 if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
843 if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
845 // ************ Eigenvector relative to value **************
847 x2 = (A33*A21-A23*A31-value*A21)/(A22*A31-A32*A21-value*A31);
848 x1 = (value-A33-A32*x2)/A31;
849 Vec.SetXYZ(x1,x2,x3);
854 nr = -n1*xbar-n2*ybar-n3*zbar;
856 u0 = -0.5*n1/(nr+n3);
857 v0 = -0.5*n2/(nr+n3);
858 rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
861 fFit += 10.*TMath::Abs(sqrt((original[0]->X()-u0)*(original[0]->X()-u0)+(original[0]->Y()-v0)*(original[0]->Y()-v0))-rho);
862 fFit += 10.*TMath::Abs(sqrt((original[1]->X()-u0)*(original[1]->X()-u0)+(original[1]->Y()-v0)*(original[1]->Y()-v0))-rho);
863 fFit += 10.*TMath::Abs(sqrt((original[2]->X()-u0)*(original[2]->X()-u0)+(original[2]->Y()-v0)*(original[2]->Y()-v0))-rho);
864 fFit += 10.*TMath::Abs(sqrt((original[3]->X()-u0)*(original[3]->X()-u0)+(original[3]->Y()-v0)*(original[3]->Y()-v0))-rho);
865 fFit += 10.*TMath::Abs(sqrt((original[4]->X()-u0)*(original[4]->X()-u0)+(original[4]->Y()-v0)*(original[4]->Y()-v0))-rho);
866 fFit += 10.*TMath::Abs(sqrt((original[5]->X()-u0)*(original[5]->X()-u0)+(original[5]->Y()-v0)*(original[5]->Y()-v0))-rho);
868 fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
869 fphi = TMath::ATan2(v0,u0);
871 //**************************************************************************
872 // LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
873 // strictly linear (no approximation)
874 //**************************************************************************
876 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
878 // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
879 // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
881 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
886 ierrl=FitLinear(N,original,linerrors,omega,u0,v0,fphi,zData,zError,CorrLin);
888 // fprintf(pout,"%f \n",chisql);
892 ierr = (ierrl > ierr ? ierrl : ierr);
898 //-----------------------------------------------------------------------------
899 void AliITSRiemannFit::Streamer(TBuffer &lRb){
900 ////////////////////////////////////////////////////////////////////////
901 // The default Streamer function "written by ROOT" doesn't write out
902 // the arrays referenced by pointers. Therefore, a specific Streamer function
903 // has to be written. This function should not be modified but instead added
904 // on to so that older versions can still be read.
905 ////////////////////////////////////////////////////////////////////////
906 // Stream an object of class AliITSRiemannFit.
911 if (lRb.IsReading()) {
912 Version_t lRv = lRb.ReadVersion(); if (lRv) { }
913 TObject::Streamer(lRb);
915 lRb >> fPrimaryTracks;
917 for(i=0;i<6;i++) lRb >> fPLay[i];
919 for(i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
921 } // end if fPointRecs!=0
922 fPointRecs = new Point_tl*[fSizeEvent];
923 for(i=0;i<fSizeEvent;i++){
924 fPointRecs[i] = new Point_tl[n];
926 lRb >> fPointRecs[i][j].lay;
927 lRb >> fPointRecs[i][j].lad;
928 lRb >> fPointRecs[i][j].det;
929 lRb >> fPointRecs[i][j].track;
930 lRb >> fPointRecs[i][j].fx;
931 lRb >> fPointRecs[i][j].fy;
932 lRb >> fPointRecs[i][j].fz;
933 lRb >> fPointRecs[i][j].fr;
934 lRb >> fPointRecs[i][j].fdE;
935 lRb >> fPointRecs[i][j].fdx;
936 lRb >> fPointRecs[i][j].fdy;
937 lRb >> fPointRecs[i][j].fdz;
938 lRb >> fPointRecs[i][j].fPt;
939 lRb >> (Char_t*)fPointRecs[i][j].fName;
941 lRb << fPointRecs[i][j].fOrigin[ii];
943 lRb << fPointRecs[i][j].fMomentum[jj];
944 lRb >> fPointRecs[i][j].fCode;
945 lRb >> fPointRecs[i][j].phi;
946 lRb >> fPointRecs[i][j].eta;
947 lRb >> fPointRecs[i][j].vertexPhi;
951 // for(i=0;i<fSizeEvent/6;i++) delete[] fspdi[i];
953 // } // end if fspdi!=0
954 // fspdi = new Point_tl*[fSizeEvent/6];
955 // for(i=0;i<fSizeEvent/6;i++){
956 // fspdi[i] = new Point_tl[n];
958 // lRb >> fspdi[i][j].lay;
959 // lRb >> fspdi[i][j].lad;
960 // lRb >> fspdi[i][j].det;
961 // lRb >> fspdi[i][j].track;
962 // lRb >> fspdi[i][j].fx;
963 // lRb >> fspdi[i][j].fy;
964 // lRb >> fspdi[i][j].fz;
965 // lRb >> fspdi[i][j].fr;
966 // lRb >> fspdi[i][j].fdE;
967 // lRb >> fspdi[i][j].fdx;
968 // lRb >> fspdi[i][j].fdy;
969 // lRb >> fspdi[i][j].fdz;
970 // lRb >> fspdi[i][j].fPt;
971 // for (ii=0;ii<4;ii++)
972 // lRb << fspdi[i][j].fOrigin[ii];
973 // for (jj=0;jj<4;jj++)
974 // lRb << fspdi[i][j].fMomentum[jj];
975 // lRb >> fspdi[i][j].fCode;
976 // lRb >> (Char_t*)fspdi[i][j].fName;
977 // lRb >> fspdi[i][j].phi;
978 // lRb >> fspdi[i][j].eta;
979 // lRb >> fspdi[i][j].vertexPhi;
983 // for(i=0;i<fSizeEvent/6;i++) delete[] fspdo[i];
985 // } // end if fspdo!=0
986 // fspdo = new Point_tl*[fSizeEvent/6];
987 // for(i=0;i<fSizeEvent/6;i++){
988 // fspdo[i] = new Point_tl[n];
990 // lRb >> fspdo[i][j].lay;
991 // lRb >> fspdo[i][j].lad;
992 // lRb >> fspdo[i][j].det;
993 // lRb >> fspdo[i][j].track;
994 // lRb >> fspdo[i][j].fx;
995 // lRb >> fspdo[i][j].fy;
996 // lRb >> fspdo[i][j].fz;
997 // lRb >> fspdo[i][j].fr;
998 // lRb >> fspdo[i][j].fdE;
999 // lRb >> fspdo[i][j].fdx;
1000 // lRb >> fspdo[i][j].fdy;
1001 // lRb >> fspdo[i][j].fdz;
1002 // lRb >> fspdo[i][j].fPt;
1003 // for (ii=0;ii<4;ii++)
1004 // lRb << fspdo[i][j].fOrigin[ii];
1005 // for (jj=0;jj<4;jj++)
1006 // lRb << fspdo[i][j].fMomentum[jj];
1007 // lRb >> fspdo[i][j].fCode;
1008 // lRb >> (Char_t*)fspdo[i][j].fName;
1009 // lRb >> fspdo[i][j].phi;
1010 // lRb >> fspdo[i][j].eta;
1011 // lRb >> fspdo[i][j].vertexPhi;
1015 lRb.WriteVersion(AliITSRiemannFit::IsA());
1016 TObject::Streamer(lRb);
1018 lRb << fPrimaryTracks;
1020 for(i=0;i<6;i++) lRb >> fPLay[i];
1021 for(i=0;i<fSizeEvent;i++) for(j=0;j<n;j++){
1022 lRb << fPointRecs[i][j].lay;
1023 lRb << fPointRecs[i][j].lad;
1024 lRb << fPointRecs[i][j].det;
1025 lRb << fPointRecs[i][j].track;
1026 lRb << fPointRecs[i][j].fx;
1027 lRb << fPointRecs[i][j].fy;
1028 lRb << fPointRecs[i][j].fz;
1029 lRb << fPointRecs[i][j].fr;
1030 lRb << fPointRecs[i][j].fdE;
1031 lRb << fPointRecs[i][j].fdx;
1032 lRb << fPointRecs[i][j].fdy;
1033 lRb << fPointRecs[i][j].fdz;
1034 lRb << fPointRecs[i][j].fPt;
1035 for (ii=0;ii<4;ii++)
1036 lRb << fPointRecs[i][j].fOrigin[ii];
1037 for (jj=0;jj<4;jj++)
1038 lRb << fPointRecs[i][j].fMomentum[jj];
1039 lRb << fPointRecs[i][j].fCode;
1040 lRb << fPointRecs[i][j].fName;
1041 lRb << fPointRecs[i][j].phi;
1042 lRb << fPointRecs[i][j].eta;
1043 lRb << fPointRecs[i][j].vertexPhi;
1045 // for(i=0;i<fSizeEvent/6;i++) for(j=0;j<n;j++){
1046 // lRb << fspdi[i][j].lay;
1047 // lRb << fspdi[i][j].lad;
1048 // lRb << fspdi[i][j].det;
1049 // lRb << fspdi[i][j].track;
1050 // lRb << fspdi[i][j].fx;
1051 // lRb << fspdi[i][j].fy;
1052 // lRb << fspdi[i][j].fz;
1053 // lRb << fspdi[i][j].fr;
1054 // lRb << fspdi[i][j].fdE;
1055 // lRb << fspdi[i][j].fdx;
1056 // lRb << fspdi[i][j].fdy;
1057 // lRb << fspdi[i][j].fdz;
1058 // lRb << fspdi[i][j].fPt;
1059 // for (ii=0;ii<4;ii++)
1060 // lRb << fspdi[i][j].fOrigin[ii];
1061 // for (jj=0;jj<4;jj++)
1062 // lRb << fspdi[i][j].fMomentum[jj];
1063 // lRb << fspdi[i][j].fCode;
1064 // lRb << fspdi[i][j].fName;
1065 // lRb << fspdi[i][j].phi;
1066 // lRb << fspdi[i][j].eta;
1067 // lRb << fspdi[i][j].vertexPhi;
1069 // for(i=0;i<fSizeEvent/6;i++) for(j=0;j<n;j++){
1070 // lRb << fspdo[i][j].lay;
1071 // lRb << fspdo[i][j].lad;
1072 // lRb << fspdo[i][j].det;
1073 // lRb << fspdo[i][j].track;
1074 // lRb << fspdo[i][j].fx;
1075 // lRb << fspdo[i][j].fy;
1076 // lRb << fspdo[i][j].fz;
1077 // lRb << fspdo[i][j].fr;
1078 // lRb << fspdo[i][j].fdE;
1079 // lRb << fspdo[i][j].fdx;
1080 // lRb << fspdo[i][j].fdy;
1081 // lRb << fspdo[i][j].fdz;
1082 // lRb << fspdo[i][j].fPt;
1083 // for (ii=0;ii<4;ii++)
1084 // lRb << fspdo[i][j].fOrigin[ii];
1085 // for (jj=0;jj<4;jj++)
1086 // lRb << fspdo[i][j].fMomentum[jj];
1087 // lRb << fspdo[i][j].fCode;
1088 // lRb << fspdo[i][j].fName;
1089 // lRb << fspdo[i][j].phi;
1090 // lRb << fspdo[i][j].eta;
1091 // lRb << fspdo[i][j].vertexPhi;