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 **************************************************************************/
20 // This class performs a fast fit of helices going through the <=6 *
21 // points of the ITS, with the goal of studying tracking and *
22 // vertexing performances. *
23 // Generated kinematics is used to take into account different weights *
24 // associated to points in different layers (with different multiple *
25 // scattering-originated errors). *
27 // Based on the work by A. Strandlie, R. Fruhwirth *
29 // First implementation by N. Bustreo, R. Turrisi - July 2000 *
31 // Further modifications by A. Dainese, R. Turrisi *
33 // Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it *
35 // ************************************************************************
38 // Modified November, 7th 2001 by Rosario Turrisi
39 // (rosario.turrisi@pd.infn.it)
41 // FitHelix returns different values. 0=ok, >0 =problem
42 // void FitLinear -> Int_t FitLinear to give feedback of errors to FitHelix
45 // Modified July, 30th 2001 by Rosario Turrisi
46 // (rosario.turrisi@pd.infn.it)
48 // Fit for z now in (z,s) plane.
49 // Returns parameters in order to write the helix equation
50 // and find the right phase/initial point.
52 // "PROPER WEIGHTS": (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2)
57 #include "AliITSRiemannFit.h"
59 #include "TClonesArray.h"
62 #include "Riostream.h"
64 #include "TGraphErrors.h"
65 #include "TParticle.h"
68 #include "AliITSRecPoint.h"
69 #include "AliITSgeom.h"
71 #include "AliITSDetTypeRec.h"
74 ClassImp(AliITSRiemannFit)
77 AliITSRiemannFit::AliITSRiemannFit():
83 ///////////////////////////////////////////////////////////
84 // Default constructor.
85 // Set everything to zero.
86 ////////////////////////////////////////////////////////////
89 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
93 //______________________________________________________________________
94 AliITSRiemannFit::AliITSRiemannFit(const AliITSRiemannFit &rf) : TObject(rf),
95 fSizeEvent(rf.fSizeEvent),
96 fPrimaryTracks(rf.fPrimaryTracks),
98 fParticles(rf.fParticles),
99 fPointRecs(rf.fPointRecs) {
104 //______________________________________________________________________
105 AliITSRiemannFit& AliITSRiemannFit::operator=(const AliITSRiemannFit& rf){
106 // Assignment operator
107 this->~AliITSRiemannFit();
108 new(this) AliITSRiemannFit(rf);
113 //______________________________________________________________________
114 AliITSRiemannFit::~AliITSRiemannFit() {
115 ///////////////////////////////////////////////////////////
116 // Default destructor.
117 // if arrays exist delete them. Then set everything to zero.
118 ////////////////////////////////////////////////////////////
120 for(Int_t i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
122 } // end if fPointRecs!=0
131 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
134 //----------------------------------------------------------------------
136 AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks):
138 fPrimaryTracks(ntracks),
142 ///////////////////////////////////////////////////////////
144 // Set fSizeEvent to size and fPrimaryTracks to ntracks.
146 ////////////////////////////////////////////////////////////
149 AliPointtl *first = new AliPointtl[fSizeEvent];
150 AliPointtl **pointRecs = new AliPointtl*[fSizeEvent];
151 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
152 for(Int_t j=0;j<fSizeEvent;j++) // create an array of struct
153 pointRecs[j] = &(first[j]);
156 // ---------------------------------------------------------------------
157 AliITSRiemannFit::AliPointtl::AliPointtl():
176 fEta(0),fVertexPhi(0){
177 // default constructor
200 AliITSRiemannFit::AliPointtl::AliPointtl(const AliPointtl& ap):
214 fMomentum(ap.fMomentum),
220 fVertexPhi(ap.fVertexPhi){
227 // ---------------------------------------------------------------------
229 void FillPoints(AliITSRiemannFit::AliPointtl **Points,Int_t &index,Float_t *xpoint,
231 TLorentzVector pE,TLorentzVector oT,Int_t *id,
232 Int_t track, Char_t *name,Int_t code,
234 ///////////////////////////////////////////////////////////////////////
235 // Fill the structure AliPointtl with the proper data
237 //////////////////////////////////////////////////////////////////////
238 Float_t pPI2 = 2.0*TMath::Pi();
246 phi = TMath::ATan2(y,x);
247 if(phi<0.0) phi += pPI2;
248 Points[i]->SetPhi(phi);
249 Points[i]->SetEta(-0.5*tan(0.5*TMath::ATan2(r,z)));
253 Points[i]->SetdX(error[0]);
254 Points[i]->SetdY(error[1]);
255 Points[i]->SetdZ(error[2]);
257 Points[i]->SetTrack(track);
258 Points[i]->SetLay(id[0]);
259 Points[i]->SetLad(id[1]);
260 Points[i]->SetDet(id[2]);
261 Points[i]->SetMomentum(&pE);
262 Points[i]->SetOrigin(&oT);
263 Points[i]->SetPt(sqrt(pE.X()*pE.X()+pE.Y()*pE.Y()));
264 Points[i]->SetCode(code);
265 Points[i]->SetName(name);
266 Points[i]->SetVertexPhi(phiorigin);
271 // -----------------------------------------------------------------------
273 void AliITSRiemannFit::InitPoints(Int_t ntracks,TTree *TR,Int_t nparticles){
274 //////////////////////////////////////////////////////////////////////
275 // Fill the class member fPointRecs with the reconstructed points
276 // Set All other members to the real values
278 /////////////////////////////////////////////////////////////////////
279 printf("\n ************* Starting Init Points *************\n");
282 AliRunLoader* rl = AliRunLoader::Open("galice.root");
284 AliITSLoader* loader = static_cast<AliITSLoader*>(rl->GetLoader("ITSLoader"));
286 Error("InitPoints", "ITS loader not found");
289 AliITSgeom* gm = loader->GetITSgeom();
291 //get pointer to modules array
292 Int_t nmodules = gm->GetIndexMax();
293 // Get the points from points file
296 AliITSRecPoint *recp;
297 nent=TR->GetEntries();
298 AliITSDetTypeRec detTypeRec;
299 TClonesArray *iTSrec = detTypeRec.RecPoints();
301 for (mod=0; mod<nmodules; mod++) {
302 detTypeRec.ResetRecPoints();
304 Int_t nrecp = iTSrec->GetEntries();
310 fPrimaryTracks = ntracks;
311 fParticles = nparticles;
312 AliITSRiemannFit::AliPointtl *global = new AliPointtl[iMAX];
313 fPointRecs = new AliITSRiemannFit::AliPointtl*[iMAX];
315 for(Int_t j=0;j<iMAX;j++) {
316 fPointRecs[j] = &(global[j]);
319 Int_t ieta=0,ieta2=0;
320 Int_t i,id[4],idold[4];
321 Int_t track=0;// // track of hit
322 Float_t xpoint[3],errorPlus[3],errorMinus[3],globalError[3]; // position and error of the point
323 TLorentzVector oT,pE;
324 Float_t locals[3],localserror[3],localsplus[3],localsminus[3]; // local position and local errors
328 Int_t layer,ladder,detector;
329 Float_t xcluster,zcluster;
330 Int_t num=0,nspdi=0,nspdo=0,nsddi=0,nsddo=0,nssdi=0,nssdo=0;
332 for (mod=0; mod<nmodules; mod++) {
333 //itsModule=(AliITSmodule*)iTSmodules->At(mod);
334 //ITS->ResetRecPoints();
335 detTypeRec.ResetRecPoints();
337 Int_t nrecp = iTSrec->GetEntries();
338 if (!nrecp) continue;
339 //itsModule->GetID(layer,ladder,detector);
340 gm->GetModuleId(mod,layer,ladder,detector);
342 for (irec=0;irec<nrecp;irec++) {
343 recp = (AliITSRecPoint*)iTSrec->UncheckedAt(irec);
344 track=recp->GetLabel(0);
345 if(track <0 ) continue;
346 xcluster=recp->GetDetLocalX(); // x on cluster
347 zcluster=recp->GetDetLocalZ(); // z on cluster
348 part = (TParticle*) gAlice->GetMCApp()->Particle(track);
349 part->ProductionVertex(oT); // set the vertex
350 part->Momentum(pE); // set the vertex momentum
351 name = part->GetName();
353 sprintf(nam2,"%s",name);
354 code = part->GetPdgCode();
360 locals[0]=xcluster; // x on cluster
361 locals[1]=0.0; // y on cluster
362 locals[2]=zcluster; // z on cluster
363 localserror[0]=sqrt(recp->GetSigmaDetLocX2());
365 localserror[2]=sqrt(recp->GetSigmaZ2());
366 localsplus[0]=xcluster+sqrt(recp->GetSigmaDetLocX2()); // x on cluster
367 if(layer==1||layer==2) localsplus[1]=0.0150/2; // y on cluster
368 else if(layer==3||layer==4) localsplus[1]=0.0280/2; // y on cluster
369 else if(layer==5||layer==6) localsplus[1]=0.0300/2; // y on cluster
370 localsplus[2]=zcluster+sqrt(recp->GetSigmaZ2()); // z on cluster
371 localsminus[0]=xcluster-sqrt(recp->GetSigmaDetLocX2()); // x on cluster
372 localsminus[1]=0.0; // y on cluster
373 localsminus[2]=zcluster-sqrt(recp->GetSigmaZ2()); // z on cluster
375 gm->LtoG(layer,ladder,detector,locals,xpoint);
376 gm->LtoG(layer,ladder,detector,localsplus,errorPlus);
377 gm->LtoG(layer,ladder,detector,localsminus,errorMinus);
378 globalError[0]=0.5*TMath::Abs(errorPlus[0]-errorMinus[0]);
379 globalError[1]=0.5*TMath::Abs(errorPlus[1]-errorMinus[1]);
380 globalError[2]=0.5*TMath::Abs(errorPlus[2]-errorMinus[2]);
382 if(TMath::Abs(part->Eta())<=1.0) ieta++;
383 if(TMath::Abs(part->Eta())<=0.5) ieta2++;
385 if(!(id[0]==idold[0]&&id[1]==idold[1]&&
386 id[2]==idold[2]&&id[3]==idold[3])) {
387 FillPoints(fPointRecs,num,xpoint,globalError,pE,oT,id,track,nam2,code,pPhi);
411 // FillPoints(fspdi,nspdi,xpoint,globalError,pE,oT,id,track,name,code,pPhi);
415 // FillPoints(fspdo,nspdo,xpoint,globalError,pE,oT,id,track,name,code,pPhi);
429 for(i=0;i<4;i++) idold[i] = id[i];
430 for(i=0;i<3;i++) xpoint[i] = 0.0;
431 } // end if id != idold
445 printf("%d primary tracks in eta=+-1\n",ieta);
446 printf("%d primary tracks#2 in eta=+-0.5\n",ieta2);
447 printf("\nInitPoints :\n\nPoints on Layer1 : %d on Layer2 : %d\n",nspdi,nspdo);
448 printf("Points on Layer3 : %d on Layer4 : %d\n",nsddi,nsddo);
449 printf("Points on Layer5 : %d on Layer6 : %d\n",nssdi,nssdo);
450 printf("Points on all Layers: %d\n",num);
451 printf("\n ************* Init Points Finished *************\n");
454 // ------------------------------------------------------------------------
455 ///////////////////////////////////////////////////////////
456 // Functions for sorting the fPointRecs array
457 ///////////////////////////////////////////////////////////
458 Bool_t SortZ(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){
459 // Z sorting function for qsort.
462 a = s1->GetZ() - s2->GetZ();
463 if(a<0.0) return kTRUE;
464 if(a>0.0) return kFALSE;
467 Bool_t SortTrack(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){
468 // track sorting function for qsort.
471 a = s1->GetTrack() - s2->GetTrack();
472 if(a<0.0) return kTRUE;
473 if(a>0.0) return kFALSE;
476 void hpsortTrack(AliITSRiemannFit::AliPointtl **ra,Int_t n){
478 AliITSRiemannFit::AliPointtl *rra;
482 l = ((n-1) >> 1) +1; // divide 2 + 1
486 rra = ra[--l]; // decrement first
490 if(--ir == 0){ // decrement first
498 if( j<ir && SortTrack(ra[j],ra[j+1]) ) j++;
499 if( SortTrack(rra,ra[j]) ){
510 void hpsortZ(AliITSRiemannFit::AliPointtl **ra,Int_t n){
512 AliITSRiemannFit::AliPointtl *rra;
516 l = ((n-1) >> 1) +1; // devide 2 + 1
520 rra = ra[--l]; // decrament first
524 if(--ir == 0){ // decrament first
532 if( j<ir && SortZ(ra[j],ra[j+1]) ) j++;
533 if( SortZ(rra,ra[j]) ){
544 //-----------------------------------------------------------------------
545 ////////////////////////////////////////////////////////////////////
547 ///////////////////////////////////////////////////////////////////
548 Int_t Partition(Int_t array[],Int_t left,Int_t right){
549 Int_t val = array[left];
551 Int_t rm = right + 1;
562 Int_t tempr = array[rm];
573 ///////////////////////////////////////////////////////////////////////
575 void AliITSRiemannFit::WritePoints(void) {
576 /////////////////////////////////////////////////////////////////////
577 // write the data in a file (temporary ascii)
578 /////////////////////////////////////////////////////////////////////
579 FILE *ascii= fopen("AsciiPoints.dat","w");
580 for(Int_t i=0;i<fPoints;i++) {
581 fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->GetLay(),
582 fPointRecs[i]->GetTrack(),fPointRecs[i]->GetX(),
583 fPointRecs[i]->GetY(),fPointRecs[i]->GetZ());
588 //-----------------------------------------------------------------------
590 void AliITSRiemannFit::ReadPoints(void) {
591 //////////////////////////////////////////////////////////////////////
592 // read the filled array
593 /////////////////////////////////////////////////////////////////////
594 hpsortTrack(fPointRecs,fPoints);
595 for(Int_t i=0;i<fPoints;i++)
596 printf("%d\t%d\t%d\t%f\t%f\t%f\t(%.0f,%.0f,%.0f)\t%.3f\t%s\n",
597 i,fPointRecs[i]->GetLay(),fPointRecs[i]->GetTrack(),
598 fPointRecs[i]->GetX(),fPointRecs[i]->GetY(),
599 fPointRecs[i]->GetZ(),fPointRecs[i]->GetOrigin()->X(),
600 fPointRecs[i]->GetOrigin()->Y(),fPointRecs[i]->GetOrigin()->Z(),
601 fPointRecs[i]->GetPt(),fPointRecs[i]->GetName());
604 //-----------------------------------------------------------------------
606 Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c,
607 Double_t &x1,Double_t &x2,Double_t &x3){
608 //////////////////////////////////////////////
609 /// Solve cubic equation:
610 /// x^3 + a*x^2 +b*x + c
612 /// returns x1 , x2 , x3
613 ////////////////////////////////////////
615 Double_t qQ = ((a*a - 3*b)/9);
616 Double_t rR = ((2*a*a*a - 9*a*b +27*c)/54);
618 Double_t aF = -2*sqrt(qQ);
620 Double_t pPI2 = TMath::Pi()*2;
622 if( rR*rR>qQ*qQ*qQ ) {
623 cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<<endl;
630 theta = TMath::ACos(rR/sqrt(qQ*qQ*qQ));
632 x1 = (aF*TMath::Cos(theta/3))-g;
633 x2 = (aF*TMath::Cos((theta+pPI2)/3))-g;
634 x3 = (aF*TMath::Cos((theta-pPI2)/3))-g;
638 //-----------------------------------------------------------------
640 void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) {
641 ///////////////////////////////////////////////////////////////////////
642 // This function apllies the transformation in the Riemann sphere
644 ///////////////////////////////////////////////////////////////////////
645 Float_t *rR = new Float_t[npoints];
646 Float_t *theta = new Float_t[npoints];
647 Float_t pPI2 = 2*TMath::Pi();
650 for(Int_t i=0;i<npoints;i++) {
651 rR[i] = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y());
652 theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X());
653 if(theta[i]<0) theta[i]+=pPI2;
654 x = rR[i]*cos(theta[i])/(1+rR[i]*rR[i]);
655 y = rR[i]*sin(theta[i])/(1+rR[i]*rR[i]);
656 z = rR[i]*rR[i]/(1+rR[i]*rR[i]);
657 To[i]->SetXYZ(x,y,z);
665 //---------------------------------------------------------------------
667 Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t omega,
668 Double_t &thu0, Double_t &thv0, Double_t &phi, TVector2 &zData, TVector3 &zError,
670 ///////////////////////////////////////////////////////////////////////
671 // Fit the points in the (z,s) plane - helix 3rd equation
673 ///////////////////////////////////////////////////////////////////////
675 //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
676 //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
677 Double_t * z = new Double_t[npoints];
678 Double_t * x = new Double_t[npoints];
679 Double_t * y = new Double_t[npoints];
680 Double_t * s = new Double_t[npoints];
681 Double_t * ez = new Double_t[npoints];
682 Double_t * ex = new Double_t[npoints];
683 Double_t * ey = new Double_t[npoints];
684 Double_t * es = new Double_t[npoints];
685 Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
687 // Double_t chi=TMath::Pi()/2.0+phi;
688 Double_t chi=-TMath::Pi()-phi;
689 Double_t angold=0.0, tpang=0.0;
690 for(Int_t k = 0; k<npoints; k++) {
691 x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
692 y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
693 z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
694 if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
696 cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
707 Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
708 if( (x[k]-thu0)<0 ) {
710 tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
715 if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
716 s[k] = (ang1+chi)/omega;
717 es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
719 if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
721 TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez);
722 fitHist->Fit("pol1","qQ");
723 z0 = fitHist->GetFunction("pol1")->GetParameter(0);
724 vpar = fitHist->GetFunction("pol1")->GetParameter(1);
725 ez0 = fitHist->GetFunction("pol1")->GetParError(0);
726 evpar = fitHist->GetFunction("pol1")->GetParError(1);
727 chisquare = fitHist->GetFunction("pol1")->GetChisquare();
729 zError.SetXYZ(ez0,evpar,chisquare);
737 for(Int_t j = 0; j < npoints; j++) {
742 avs /= (Double_t)npoints;
743 avz /= (Double_t)npoints;
744 avsz /= (Double_t)npoints;
746 for(Int_t l = 0; l < npoints; l++) {
747 sigmas += (s[l]-avs)*(s[l]-avs);
748 sigmaz += (z[l]-avz)*(z[l]-avz);
750 sigmas /=(Double_t)npoints;
751 sigmaz /=(Double_t)npoints;
753 sigmas = sqrt(sigmas);
754 sigmaz = sqrt(sigmaz);
756 corrLin = (avsz-avs*avz)/(sigmas*sigmaz);
770 //-------------------------------------------------------------------
771 Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
772 Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,Double_t& omega, Double_t& z0,
773 Double_t& vpar,Double_t& chisql, Double_t& fCorrLin,Double_t& fFit,
774 Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
775 ///////////////////////////////////////////////////////////////////////
776 // This function finds the helix paramenters
777 // d0 = impact parameter
778 // rho = radius of circle
781 // starting from the momentum and the outcome of
782 // the fit on the Riemann sphere (i.e. u0,v0,rho)
784 // MIND !!!! Here we assume both angular velocities be 1.0 (yes, one-dot-zero !)
787 ///////////////////////////////////////////////////////////////////////
789 // All this stuff relies on this hypothesis !!!
791 // FILE *pout=fopen("chisql.dat","a");
792 Int_t ierr = 0, ierrl=0;
795 Int_t bitlay[6]={1,1,1,1,1,1};
796 bitlay[0]*=first; bitlay[1]*=second; bitlay[2]*=third; bitlay[3]*=fourth; bitlay[4]*=fifth; bitlay[5]*=sixth;
797 fd0 = -9999; // No phisycs value
798 u0 = -9999.9999; // parameters of helix - strange value...
799 v0 = -9999.9999; // parameters of helix - strange value...
800 rho = -9999.9999; // parameters of helix -unphysical strange value...
802 const Char_t* name = 0;
806 Int_t npl[6]={0,0,0,0,0,0};
807 Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
808 Double_t pt = sqrt(Px*Px+Py*Py);
812 TVector3 *ori = new TVector3[iMAX];
813 TVector3 **original = new TVector3*[iMAX];
814 TVector3 *rie = new TVector3[iMAX];
815 TVector3 **riemann = new TVector3*[iMAX];
816 TVector3 *err = new TVector3[iMAX];
817 TVector3 **errors = new TVector3*[iMAX];
818 TVector3 *linerr = new TVector3[iMAX];
819 TVector3 **linerrors = new TVector3*[iMAX];
820 //PH Double_t weight[iMAX];
821 Double_t * weight = new Double_t[iMAX];
824 original[i] = &(ori[i]);
825 riemann[i] = &(rie[i]);
826 errors[i] = &(err[i]);
827 linerrors[i] = &(linerr[i]);
829 for(k =0;k<iMAX;k++) original[k]->SetXYZ(9999,9999,9999);
830 Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
831 a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
835 Double_t a,b,c,d; // cubic parameters
836 Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
837 Double_t value = 0.0; // minimum eigenvalue
838 Double_t x1,x2,x3; // eigenvector component
839 Double_t n1,n2,n3,nr= 0;// unit eigenvector
840 Double_t radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm]
841 Double_t sigmaMS = 0;
842 TVector3 vVec,vVecNor;
844 // Select RecPoints belonging to the track
845 for(k =0;k<fPoints;k++){
846 if(fPointRecs[k]->GetTrack()==tracknumber) {
847 name = fPointRecs[k]->GetName();
848 pt = fPointRecs[k]->GetPt();
849 pLayer = fPointRecs[k]->GetLay();
850 Int_t ilay = pLayer-1;
851 if(npl[ilay]!=0) continue;
852 if(bitlay[ilay] == 1) {
853 original[nN]->SetXYZ(0.1*fPointRecs[k]->GetX(),0.1*fPointRecs[k]->GetY(),0.1*fPointRecs[k]->GetZ());
854 errors[nN]->SetXYZ(0.1*fPointRecs[k]->GetdX(),0.1*fPointRecs[k]->GetdY(),0.1*fPointRecs[k]->GetdZ());
855 sigmaMS = (radiusdm[pLayer]-radiusdm[0])*0.000724/pP;// beam pipe contribution
856 for(Int_t j=1;j<pLayer;j++) {
857 sigmaMS += (radiusdm[pLayer]-radiusdm[j])*0.00136/pP;
859 weight[nN] = ( 1 + original[nN]->Perp2() )*( 1+ original[nN]->Perp2() )/
860 ( errors[nN]->Perp2() + sigmaMS*sigmaMS );
861 linerrors[nN]->SetXYZ(errors[nN]->X(),errors[nN]->Y(),sqrt(errors[nN]->Z()*errors[nN]->Z()+sigmaMS*sigmaMS));
865 } //end if track==tracknumber
868 // 6 points, no more, no less
870 if(original[5]->X() == 9999 || original[6]->X() != 9999)
873 return 1; // not enough points
879 // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
882 RiemannTransf(nN,original,riemann);
884 Double_t sumWeights = 0.0; // sum of weights factor
886 for(Int_t j=0;j<nN;j++){ // mean values for x[i],y[i],z[i]
887 xbar+=weight[j]*riemann[j]->X();
888 ybar+=weight[j]*riemann[j]->Y();
889 zbar+=weight[j]*riemann[j]->Z();
890 sumWeights+=weight[j];
897 for(Int_t j=0;j<nN;j++) { // Calculate the matrix elements
898 a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
899 a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
900 a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
901 a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
902 a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
903 a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
916 // ************** Determinant parameters ********************
917 // n.b. simplifications done keeping in mind symmetry of A
921 c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
922 d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
924 // ************** Find the 3 eigenvalues *************************
925 Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
928 printf("Track %d Has no real solution continuing ...\n",tracknumber);
933 // **************** Find the lowest eigenvalue *****************
934 if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
935 if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
936 if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
938 // ************ Eigenvector relative to value **************
940 x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
941 x1 = (value-a33-a32*x2)/a31;
942 vVec.SetXYZ(x1,x2,x3);
943 vVecNor = vVec.Unit();
947 nr = -n1*xbar-n2*ybar-n3*zbar;
949 u0 = -0.5*n1/(nr+n3);
950 v0 = -0.5*n2/(nr+n3);
951 rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
954 fFit += 10.*TMath::Abs(sqrt((original[0]->X()-u0)*(original[0]->X()-u0)+(original[0]->Y()-v0)*(original[0]->Y()-v0))-rho);
955 fFit += 10.*TMath::Abs(sqrt((original[1]->X()-u0)*(original[1]->X()-u0)+(original[1]->Y()-v0)*(original[1]->Y()-v0))-rho);
956 fFit += 10.*TMath::Abs(sqrt((original[2]->X()-u0)*(original[2]->X()-u0)+(original[2]->Y()-v0)*(original[2]->Y()-v0))-rho);
957 fFit += 10.*TMath::Abs(sqrt((original[3]->X()-u0)*(original[3]->X()-u0)+(original[3]->Y()-v0)*(original[3]->Y()-v0))-rho);
958 fFit += 10.*TMath::Abs(sqrt((original[4]->X()-u0)*(original[4]->X()-u0)+(original[4]->Y()-v0)*(original[4]->Y()-v0))-rho);
959 fFit += 10.*TMath::Abs(sqrt((original[5]->X()-u0)*(original[5]->X()-u0)+(original[5]->Y()-v0)*(original[5]->Y()-v0))-rho);
961 fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
962 fphi = TMath::ATan2(v0,u0);
964 //**************************************************************************
965 // LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
966 // strictly linear (no approximation)
967 //**************************************************************************
969 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
971 // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
972 // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
974 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
979 ierrl=FitLinear(nN,original,linerrors,omega,u0,v0,fphi,zData,zError,corrLin);
981 // fprintf(pout,"%f \n",chisql);
985 ierr = (ierrl > ierr ? ierrl : ierr);
990 Int_t AliITSRiemannFit::FitHelix(Int_t NPoints, TVector3** fPointRecs,TVector3** fPointRecErrors,Float_t& f1, Float_t& f2, Float_t& f3) {
992 ///////////////////////////////////////////////////////////////////////
993 // This function finds the helix parameters
994 // d0 = impact parameter
995 // rho = radius of circle
998 // starting from the momentum and the outcome of
999 // the fit on the Riemann sphere (i.e. u0,v0,rho)
1001 // MIND !!!! Here we assume both angular velocities be 1.0e-2 (yes, 0.01 !)
1004 // Also linear fit in (z,s) is performed, so it's 3-D !
1005 // z0 and vpar are calculated (intercept and z-component of velocity, but
1006 // in units... you guess.
1009 // Values calculated in addition:
1011 // - transverse impact parameter fd0
1012 // - sum of residuals in (x,y) plane fFit
1013 // - chisquare of linear fit chisql
1014 // - correlation coefficient fCorrLin
1020 ///////////////////////////////////////////////////////////////////////
1022 // All this stuff relies on this hypothesis !!!
1024 Int_t ierr = 0, ierrl=0;
1025 const Double_t kOmega = 1.0e-2;
1030 Double_t fd0 = -9999; // fake values
1031 Double_t u0 = -9999.9999; // for eventual
1032 Double_t v0 = -9999.9999; // debugging
1033 Double_t rho = -9999.9999; //
1034 Double_t fphi, fFit, chisql, z0, vpar, fCorrLin;
1037 // This info is no more there... to be re-considered... maybe
1039 // Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
1040 // Double_t pt = sqrt(Px*Px+Py*Py);
1045 TVector3 *ori = new TVector3[NPoints];
1046 TVector3 **original = new TVector3*[NPoints];
1047 TVector3 *rie = new TVector3[NPoints];
1048 TVector3 **riemann = new TVector3*[NPoints];
1049 TVector3 *err = new TVector3[NPoints];
1050 TVector3 **errors = new TVector3*[NPoints];
1051 TVector3 *linerr = new TVector3[NPoints];
1052 TVector3 **linerrors = new TVector3*[NPoints];
1053 Double_t * weight = new Double_t[NPoints];
1055 for(Int_t i=0; i<NPoints; i++){
1057 original[i] = &(ori[i]);
1058 riemann[i] = &(rie[i]);
1059 errors[i] = &(err[i]);
1060 linerrors[i] = &(linerr[i]);
1062 original[i]->SetXYZ(9999,9999,9999);
1066 // Riemann fit parameters
1068 Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
1069 a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
1074 Double_t a,b,c,d; // cubic parameters
1075 Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
1076 Double_t value = 0.0; // minimum eigenvalue
1077 Double_t x1,x2,x3; // eigenvector component
1078 Double_t n1,n2,n3,nr= 0; // unit eigenvector
1079 TVector3 vVec,vVecNor;
1081 for (Int_t ip=0; ip<NPoints; ip++) {
1082 original[ip]->SetXYZ(0.1*fPointRecs[ip]->X(),0.1*fPointRecs[ip]->Y(),0.1*fPointRecs[ip]->Z());
1084 errors[ip]->SetXYZ(0.1*fPointRecErrors[ip]->X(),0.1*fPointRecErrors[ip]->Y(),0.1*fPointRecErrors[ip]->Z());
1085 weight[ip] = (1+original[ip]->Perp2())*(1+original[ip]->Perp2())/(errors[ip]->Perp2());
1086 linerrors[ip]->SetXYZ(errors[ip]->X(),errors[ip]->Y(),errors[ip]->Z());
1092 // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
1095 RiemannTransf(NPoints,original,riemann);
1097 Double_t sumWeights = 0.0; // sum of weights factor
1099 for(Int_t j=0;j<NPoints;j++){ // mean values for x[i],y[i],z[i]
1100 xbar+=weight[j]*riemann[j]->X();
1101 ybar+=weight[j]*riemann[j]->Y();
1102 zbar+=weight[j]*riemann[j]->Z();
1103 sumWeights+=weight[j];
1110 for(Int_t j=0;j<NPoints;j++) { // Calculate the matrix elements
1111 a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
1112 a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
1113 a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
1114 a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
1115 a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
1116 a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
1119 // this doesn't seem to work...
1121 // a11 /= sumWeights;
1122 // a12 /= sumWeights;
1123 // a22 /= sumWeights;
1124 // a23 /= sumWeights;
1125 // a13 /= sumWeights;
1126 // a33 /= sumWeights;
1138 // ************** Determinant parameters ********************
1139 // n.b. simplifications done keeping in mind symmetry of A
1143 c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
1144 d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
1146 // ************** Find the 3 eigenvalues *************************
1147 Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
1149 if(checkCubic !=1 ){
1150 printf("No real solution. Check data.\n");
1155 // **************** Find the lowest eigenvalue *****************
1156 if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
1157 if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
1158 if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
1160 // ************ Eigenvector relative to value **************
1162 x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
1163 x1 = (value-a33-a32*x2)/a31;
1164 vVec.SetXYZ(x1,x2,x3);
1165 vVecNor = vVec.Unit();
1169 nr = -n1*xbar-n2*ybar-n3*zbar;
1171 u0 = -0.5*n1/(nr+n3);
1172 v0 = -0.5*n2/(nr+n3);
1173 rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
1177 for (Int_t i=0; i<NPoints; i++) {
1178 fFit += 10.*TMath::Abs(sqrt((original[i]->X()-u0)*(original[i]->X()-u0)+(original[i]->Y()-v0)*(original[i]->Y()-v0))-rho);
1180 fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
1181 fphi = TMath::ATan2(v0,u0);
1183 //**************************************************************************
1184 // LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
1185 // strictly linear (no approximation)
1186 //**************************************************************************
1188 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1190 // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
1191 // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
1193 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1199 ierrl=LinearFit(NPoints,original,linerrors,kOmega,u0,v0,fphi,zData,zError,corrLin);
1200 if(ierrl==33) return 0;
1202 // fprintf(pout,"%f \n",chisql);
1206 ierr = (ierrl > ierr ? ierrl : ierr);
1211 f2=vpar/(kOmega*TMath::Abs(rho));
1227 //____________________________________________________________
1229 Int_t AliITSRiemannFit::LinearFit(Int_t npoints, TVector3 **input,
1230 TVector3 **errors, Double_t omega,
1231 Double_t &thu0, Double_t &thv0, Double_t &phi,TVector2 &zData, TVector3 &zError,
1233 ///////////////////////////////////////////////////////////////////////
1234 // Fit the points in the (z,s) plane - helix 3rd equation
1236 ///////////////////////////////////////////////////////////////////////
1240 //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
1241 //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
1242 Double_t * z = new Double_t[npoints];
1243 Double_t * x = new Double_t[npoints];
1244 Double_t * y = new Double_t[npoints];
1245 Double_t * s = new Double_t[npoints];
1246 Double_t * ez = new Double_t[npoints];
1247 Double_t * ex = new Double_t[npoints];
1248 Double_t * ey = new Double_t[npoints];
1249 Double_t * es = new Double_t[npoints];
1250 Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
1253 // Double_t chi=TMath::Pi()/2.0+phi;
1254 Double_t chi=-TMath::Pi()-phi;
1255 Double_t angold=0.0, tpang=0.0;
1256 for(Int_t k = 0; k<npoints; k++) {
1257 x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
1258 y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
1259 z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
1260 if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
1262 cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
1273 Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
1274 if( (x[k]-thu0)<0 ) {
1275 if (ang1*angold<0) {
1276 tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
1281 if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
1282 s[k] = (ang1+chi)/omega;
1283 es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
1285 if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
1287 // if(s[0]>-636 && s[0]<-625) return 33;
1289 TGraph* fitHist = new TGraph(npoints,s,z);
1290 TF1* f1 = new TF1("f1",Fitfunction,-100,100,2);
1292 f1->SetParameter(0,1);
1293 f1->SetParameter(1,1);
1294 f1->SetLineColor(2);
1295 fitHist->Fit(f1,"qQ");
1297 z0 = f1->GetParameter(0);
1298 vpar = f1->GetParameter(1);
1299 ez0 = f1->GetParError(0);
1300 evpar= f1->GetParError(1);
1301 chisquare=f1->GetChisquare();
1303 zError.SetXYZ(ez0,evpar,chisquare);
1311 for(Int_t j = 0; j < npoints; j++) {
1316 avs /= (Double_t)npoints;
1317 avz /= (Double_t)npoints;
1318 avsz /= (Double_t)npoints;
1320 for(Int_t l = 0; l < npoints; l++) {
1321 sigmas += (s[l]-avs)*(s[l]-avs);
1322 sigmaz += (z[l]-avz)*(z[l]-avz);
1324 sigmas /=(Double_t)npoints;
1325 sigmaz /=(Double_t)npoints;
1327 sigmas = sqrt(sigmas);
1328 sigmaz = sqrt(sigmaz);
1330 corrLin = (avsz-avs*avz)/(sigmas*sigmaz);
1342 delete f1; delete fitHist;
1347 //_______________________________________________________
1349 Double_t AliITSRiemannFit::Fitfunction(Double_t *x, Double_t* par){
1350 // function used for fit
1351 return par[0]+(*x)*par[1];