]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSRiemannFit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Codinf conventions
[u/mrichter/AliRoot.git] / ITS / AliITSRiemannFit.cxx
CommitLineData
8db76038 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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. *
14 * *
15 * *
16 * ///////////////////////////////////////////////////////////////////// *
17 * *
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). *
24 * *
25 * Based on the work by A. Strandlie, R. Fruhwirth *
26 * *
27 * First implementation by N. Bustreo, R. Turrisi - July 2000 *
28 * *
29 * Further modifications by A. Dainese, R. Turrisi *
30 * *
31 * Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it *
32 * *
33 * **************************************************************************/
34//
35//
36// Modified November, 7th 2001 by Rosario Turrisi
37// (rosario.turrisi@pd.infn.it)
38//
39// FitHelix returns different values. 0=ok, >0 =problem
40// void FitLinear -> Int_t FitLinear to give feedback of errors to FitHelix
41//
42//
43// Modified July, 30th 2001 by Rosario Turrisi
44// (rosario.turrisi@pd.infn.it)
45//
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.
49//
50// "PROPER WEIGHTS": (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2)
51//
4ae5bbc4 52#include <Riostream.h>
5ff6c5cc 53#include "AliITS.h"
8db76038 54#include "AliITSRiemannFit.h"
55#include "AliRun.h"
56#include "TClonesArray.h"
57#include "stdio.h"
58#include "stdlib.h"
4ae5bbc4 59#include "Riostream.h"
8db76038 60#include "TMath.h"
61#include "TF1.h"
62#include "TGraphErrors.h"
8db76038 63#include "TStyle.h"
8db76038 64#include "TParticle.h"
5ff6c5cc 65#include "TTree.h"
66#include "TVector3.h"
8db76038 67#include "AliITSRecPoint.h"
68#include "AliITSgeom.h"
69#include "AliITSmodule.h"
5d12ce38 70#include "AliMC.h"
8db76038 71ClassImp(AliITSRiemannFit)
72
73
74AliITSRiemannFit::AliITSRiemannFit() {
75 ///////////////////////////////////////////////////////////
76 // Default constructor.
77 // Set everything to zero.
78 ////////////////////////////////////////////////////////////
79
80 fSizeEvent = 0;
81 fPoints = 0;
82 fPrimaryTracks = 0;
83 fPointRecs = 0;
84 //
85 // test erase
86// fspdi = 0;
87// fspdo = 0;
88 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
89
90}
91//----------------------------------------------------------------------
92
93AliITSRiemannFit::~AliITSRiemannFit() {
94 ///////////////////////////////////////////////////////////
95 // Default destructor.
96 // if arrays exist delete them. Then set everything to zero.
97 ////////////////////////////////////////////////////////////
98 if(fPointRecs!=0){
99 for(Int_t i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
100 delete[] fPointRecs;
101 } // end if fPointRecs!=0
102 fSizeEvent = 0;
103 fPointRecs = 0;
104 fPoints = 0;
105 fPrimaryTracks = 0;
106 //
107 // test erase
108// fspdi = 0;
109// fspdo = 0;
110 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
111 return;
112}
113//----------------------------------------------------------------------
114
115AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks) {
116 ///////////////////////////////////////////////////////////
117 // Constructor.
118 // Set fSizeEvent to size and fPrimaryTracks to ntracks.
119 // Others to zero.
120 ////////////////////////////////////////////////////////////
121
122 fSizeEvent = size;
123 fPoints = 0;
124 fPrimaryTracks = ntracks;
125 //
126 // test erase
127// fspdi = 0;
128// fspdo = 0;
5ff6c5cc 129 AliPointtl *first = new AliPointtl[fSizeEvent];
130 AliPointtl **pointRecs = new AliPointtl*[fSizeEvent];
8db76038 131 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
132 for(Int_t j=0;j<fSizeEvent;j++) // create an array of struct
5ff6c5cc 133 pointRecs[j] = &(first[j]);
134}
135
136// ---------------------------------------------------------------------
137AliITSRiemannFit::AliPointtl::AliPointtl(){
138 // default constructor
139 SetLay();
140 SetLad();
141 SetDet();
142 SetTrack();
143 SetX();
144 SetY();
145 SetZ();
146 SetR();
147 SetdE();
148 SetdX();
149 SetdY();
150 SetdZ();
151 SetOrigin();
152 SetMomentum();
153 SetCode();
154 SetName();
155 SetPt();
156 SetPhi();
157 SetEta();
158 SetVertexPhi();
8db76038 159}
5ff6c5cc 160
8db76038 161// ---------------------------------------------------------------------
162
5ff6c5cc 163void FillPoints(AliITSRiemannFit::AliPointtl **Points,Int_t &index,Float_t *xpoint,
8db76038 164 Float_t *error,
5ff6c5cc 165 TLorentzVector pE,TLorentzVector oT,Int_t *id,
166 Int_t track, Char_t *name,Int_t code,
8db76038 167 Float_t phiorigin){
168 ///////////////////////////////////////////////////////////////////////
5ff6c5cc 169 // Fill the structure AliPointtl with the proper data
8db76038 170 //
171 //////////////////////////////////////////////////////////////////////
5ff6c5cc 172 Float_t pPI2 = 2.0*TMath::Pi();
8db76038 173 Float_t phi,r,x,y,z;
174 Int_t i;
175 i = index;
176 x = xpoint[0];
177 y = xpoint[1];
178 z = xpoint[2];
179 r = sqrt(x*x+y*y);
180 phi = TMath::ATan2(y,x);
5ff6c5cc 181 if(phi<0.0) phi += pPI2;
182 Points[i]->SetPhi(phi);
183 Points[i]->SetEta(-0.5*tan(0.5*TMath::ATan2(r,z)));
184 Points[i]->SetX(x);
185 Points[i]->SetY(y);
186 Points[i]->SetZ(z);
187 Points[i]->SetdX(error[0]);
188 Points[i]->SetdY(error[1]);
189 Points[i]->SetdZ(error[2]);
190 Points[i]->SetR(r);
191 Points[i]->SetTrack(track);
192 Points[i]->SetLay(id[0]);
193 Points[i]->SetLad(id[1]);
194 Points[i]->SetDet(id[2]);
195 Points[i]->SetMomentum(&pE);
196 Points[i]->SetOrigin(&oT);
197 Points[i]->SetPt(sqrt(pE.X()*pE.X()+pE.Y()*pE.Y()));
198 Points[i]->SetCode(code);
199 Points[i]->SetName(name);
200 Points[i]->SetVertexPhi(phiorigin);
8db76038 201 index++;
202 return;
203
204}
205// -----------------------------------------------------------------------
206
088e0b8d 207void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS,
8db76038 208 TTree *TR,Int_t nparticles){
209 //////////////////////////////////////////////////////////////////////
210 // Fill the class member fPointRecs with the reconstructed points
211 // Set All other members to the real values
212 //
213 /////////////////////////////////////////////////////////////////////
214 printf("\n ************* Starting Init Points *************\n");
215 TParticle *part;
216 AliITSgeom *gm = (AliITSgeom*)ITS->GetITSgeom();
217 //get pointer to modules array
5ff6c5cc 218 TObjArray *iTSmodules = ITS->GetModules();
219 Int_t nmodules=iTSmodules->GetEntriesFast();
8db76038 220 printf("nmodules = %d \n",nmodules);
221 // Get the points from points file
222 AliITSmodule *itsModule;
223 Int_t mod,irec;
224 Stat_t nent;
225 AliITSRecPoint *recp;
226 nent=TR->GetEntries();
5ff6c5cc 227 TClonesArray *iTSrec = ITS->RecPoints();
8db76038 228
5ff6c5cc 229 Int_t totRP=0;
8db76038 230 for (mod=0; mod<nmodules; mod++) {
5ff6c5cc 231 itsModule=(AliITSmodule*)iTSmodules->At(mod);
8db76038 232 ITS->ResetRecPoints();
233 TR->GetEvent(mod);
5ff6c5cc 234 Int_t nrecp = iTSrec->GetEntries();
8db76038 235 if(!nrecp) continue;
5ff6c5cc 236 totRP += nrecp;
8db76038 237 }
238
5ff6c5cc 239 Int_t iMAX = totRP;
8db76038 240 fPrimaryTracks = ntracks;
241 fParticles = nparticles;
5ff6c5cc 242 AliITSRiemannFit::AliPointtl *global = new AliPointtl[iMAX];
243 fPointRecs = new AliITSRiemannFit::AliPointtl*[iMAX];
8db76038 244 //
8db76038 245 for(Int_t j=0;j<iMAX;j++) {
246 fPointRecs[j] = &(global[j]);
8db76038 247 }
248
249 Int_t ieta=0,ieta2=0;
250 Int_t i,id[4],idold[4];
251 Int_t track=0;// // track of hit
5ff6c5cc 252 Float_t xpoint[3],errorPlus[3],errorMinus[3],globalError[3]; // position and error of the point
253 TLorentzVector oT,pE;
254 Float_t locals[3],localserror[3],localsplus[3],localsminus[3]; // local position and local errors
255 Float_t pPhi;
8db76038 256 Int_t code;
257 const char *name;
258 Int_t layer,ladder,detector;
259 Float_t xcluster,zcluster;
260 Int_t num=0,nspdi=0,nspdo=0,nsddi=0,nsddo=0,nssdi=0,nssdo=0;
261
262 for (mod=0; mod<nmodules; mod++) {
5ff6c5cc 263 itsModule=(AliITSmodule*)iTSmodules->At(mod);
8db76038 264 ITS->ResetRecPoints();
265 TR->GetEvent(mod);
5ff6c5cc 266 Int_t nrecp = iTSrec->GetEntries();
8db76038 267 if (!nrecp) continue;
268 itsModule->GetID(layer,ladder,detector);
269
270 for (irec=0;irec<nrecp;irec++) {
5ff6c5cc 271 recp = (AliITSRecPoint*)iTSrec->UncheckedAt(irec);
8db76038 272 track=recp->fTracks[0];
273 if(track <0 ) continue;
274 xcluster=recp->GetX(); // x on cluster
275 zcluster=recp->GetZ(); // z on cluster
5d12ce38 276 part = (TParticle*) gAlice->GetMCApp()->Particle(track);
5ff6c5cc 277 part->ProductionVertex(oT); // set the vertex
278 part->Momentum(pE); // set the vertex momentum
8db76038 279 name = part->GetName();
5ff6c5cc 280 Char_t nam2[50];
281 sprintf(nam2,"%s",name);
8db76038 282 code = part->GetPdgCode();
5ff6c5cc 283 pPhi = part->Phi();
8db76038 284 id[0]=layer;
285 id[1]=ladder;
286 id[2]=detector;
287 id[3]=irec;
288 locals[0]=xcluster; // x on cluster
289 locals[1]=0.0; // y on cluster
290 locals[2]=zcluster; // z on cluster
5ff6c5cc 291 localserror[0]=sqrt(recp->GetSigmaX2());
292 localserror[1]=0.0;
293 localserror[2]=sqrt(recp->GetSigmaZ2());
294 localsplus[0]=xcluster+sqrt(recp->GetSigmaX2()); // x on cluster
295 if(layer==1||layer==2) localsplus[1]=0.0150/2; // y on cluster
296 else if(layer==3||layer==4) localsplus[1]=0.0280/2; // y on cluster
297 else if(layer==5||layer==6) localsplus[1]=0.0300/2; // y on cluster
298 localsplus[2]=zcluster+sqrt(recp->GetSigmaZ2()); // z on cluster
299 localsminus[0]=xcluster-sqrt(recp->GetSigmaX2()); // x on cluster
300 localsminus[1]=0.0; // y on cluster
301 localsminus[2]=zcluster-sqrt(recp->GetSigmaZ2()); // z on cluster
8db76038 302
303 gm->LtoG(layer,ladder,detector,locals,xpoint);
5ff6c5cc 304 gm->LtoG(layer,ladder,detector,localsplus,errorPlus);
305 gm->LtoG(layer,ladder,detector,localsminus,errorMinus);
306 globalError[0]=0.5*TMath::Abs(errorPlus[0]-errorMinus[0]);
307 globalError[1]=0.5*TMath::Abs(errorPlus[1]-errorMinus[1]);
308 globalError[2]=0.5*TMath::Abs(errorPlus[2]-errorMinus[2]);
8db76038 309 if(track<ntracks) {
310 if(TMath::Abs(part->Eta())<=1.0) ieta++;
311 if(TMath::Abs(part->Eta())<=0.5) ieta2++;
312 }
313 if(!(id[0]==idold[0]&&id[1]==idold[1]&&
314 id[2]==idold[2]&&id[3]==idold[3])) {
5ff6c5cc 315 FillPoints(fPointRecs,num,xpoint,globalError,pE,oT,id,track,nam2,code,pPhi);
8db76038 316 //
317 // test erase
318 switch (idold[0]) {
319 case 1:
320 nspdi++;
321 break;
322 case 2:
323 nspdo++;
324 break;
325 case 3:
326 nsddi++;
327 break;
328 case 4:
329 nsddo++;
330 break;
331 case 5:
332 nssdi++;
333 break;
334 case 6:
335 nssdo++;
336 break;
337 }
338// if(idold[0]==1){
5ff6c5cc 339// FillPoints(fspdi,nspdi,xpoint,globalError,pE,oT,id,track,name,code,pPhi);
8db76038 340// }
341// if(idold[0]==2){
342
5ff6c5cc 343// FillPoints(fspdo,nspdo,xpoint,globalError,pE,oT,id,track,name,code,pPhi);
8db76038 344// }
345// if(idold[0]==3){
346// nsddi++;
347// }
348// if(idold[0]==4){
349// nsddo++;
350// }
351// if(idold[0]==5){
352// nssdi++;
353// }
354// if(idold[0]==6){
355// nssdo++;
356// }
357 for(i=0;i<4;i++) idold[i] = id[i];
358 for(i=0;i<3;i++) xpoint[i] = 0.0;
359 } // end if id != idold
360 } // end for irec
361 }// end for mod
362
363 fPoints = num;
364 fSizeEvent = num;
365 fPLay[0] = nspdi ;
366 fPLay[1] = nspdo ;
367 fPLay[2] = nsddi ;
368 fPLay[3] = nsddo ;
369 fPLay[4] = nssdi ;
370 fPLay[5] = nssdo ;
371 printf("%d primary tracks in eta=+-1\n",ieta);
372 printf("%d primary tracks#2 in eta=+-0.5\n",ieta2);
373 printf("\nInitPoints :\n\nPoints on Layer1 : %d on Layer2 : %d\n",nspdi,nspdo);
374 printf("Points on Layer3 : %d on Layer4 : %d\n",nsddi,nsddo);
375 printf("Points on Layer5 : %d on Layer6 : %d\n",nssdi,nssdo);
376 printf("Points on all Layers: %d\n",num);
377 printf("\n ************* Init Points Finished *************\n");
378 return;
379}
380// ------------------------------------------------------------------------
381///////////////////////////////////////////////////////////
382// Functions for sorting the fPointRecs array
383///////////////////////////////////////////////////////////
5ff6c5cc 384Bool_t SortZ(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){
8db76038 385 // Z sorting function for qsort.
386 Float_t a;
387
5ff6c5cc 388 a = s1->GetZ() - s2->GetZ();
8db76038 389 if(a<0.0) return kTRUE;
390 if(a>0.0) return kFALSE;
391 return kFALSE;
392}
5ff6c5cc 393Bool_t SortTrack(const AliITSRiemannFit::AliPointtl *s1,const AliITSRiemannFit::AliPointtl *s2){
8db76038 394 // track sorting function for qsort.
395 Float_t a;
396
5ff6c5cc 397 a = s1->GetTrack() - s2->GetTrack();
8db76038 398 if(a<0.0) return kTRUE;
399 if(a>0.0) return kFALSE;
400 return kFALSE;
401}
5ff6c5cc 402void hpsortTrack(AliITSRiemannFit::AliPointtl **ra,Int_t n){
8db76038 403 Int_t i,ir,j,l;
5ff6c5cc 404 AliITSRiemannFit::AliPointtl *rra;
8db76038 405
406 if(n<2) return;
407
408 l = ((n-1) >> 1) +1; // divide 2 + 1
409 ir = n-1;
410 for(;;){
411 if(l>0){
412 rra = ra[--l]; // decrement first
413 }else{
414 rra = ra[ir];
415 ra[ir] = ra[0];
416 if(--ir == 0){ // decrement first
417 ra[0] = rra;
418 break;
419 } // if --ra == 0
420 } // end l>0
421 i = l;
422 j = l+1;
423 while(j<=ir){
424 if( j<ir && SortTrack(ra[j],ra[j+1]) ) j++;
425 if( SortTrack(rra,ra[j]) ){
426 ra[i] = ra[j];
427 i = j;
428 j <<= 1; // time 2.
429 }else{
430 break;
431 } // end if func()
432 } // end while
433 ra[i] = rra;
434 } // end for ever
435}
5ff6c5cc 436void hpsortZ(AliITSRiemannFit::AliPointtl **ra,Int_t n){
8db76038 437 Int_t i,ir,j,l;
5ff6c5cc 438 AliITSRiemannFit::AliPointtl *rra;
8db76038 439
440 if(n<2) return;
441
442 l = ((n-1) >> 1) +1; // devide 2 + 1
443 ir = n-1;
444 for(;;){
445 if(l>0){
446 rra = ra[--l]; // decrament first
447 }else{
448 rra = ra[ir];
449 ra[ir] = ra[0];
450 if(--ir == 0){ // decrament first
451 ra[0] = rra;
452 break;
453 } // if --ra == 0
454 } // end l>0
455 i = l;
456 j = l+1;
457 while(j<=ir){
458 if( j<ir && SortZ(ra[j],ra[j+1]) ) j++;
459 if( SortZ(rra,ra[j]) ){
460 ra[i] = ra[j];
461 i = j;
462 j <<= 1; // time 2.
463 }else{
464 break;
465 } // end if func()
466 } // end while
467 ra[i] = rra;
468 } // end for ever
469}
470//-----------------------------------------------------------------------
471////////////////////////////////////////////////////////////////////
472// Sorting functions
473///////////////////////////////////////////////////////////////////
474Int_t Partition(Int_t array[],Int_t left,Int_t right){
475 Int_t val = array[left];
476 Int_t lm = left - 1;
477 Int_t rm = right + 1;
478 for(;;) {
479 do
480 rm--;
481 while
482 (array[rm]>val);
483 do
484 lm++;
485 while
486 (array[lm]<val);
487 if(lm<rm){
488 Int_t tempr = array[rm];
489 array[rm]=array[lm];
490 array[lm]=tempr;
491 }
492 else
493 return rm;
494 }
495
496 return 1;
497}
498
499///////////////////////////////////////////////////////////////////////
500
501void AliITSRiemannFit::WritePoints(void) {
502 /////////////////////////////////////////////////////////////////////
503 // write the data in a file (temporary ascii)
504 /////////////////////////////////////////////////////////////////////
505 FILE *ascii= fopen("AsciiPoints.dat","w");
506 for(Int_t i=0;i<fPoints;i++) {
5ff6c5cc 507 fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->GetLay(),
508 fPointRecs[i]->GetTrack(),fPointRecs[i]->GetX(),
509 fPointRecs[i]->GetY(),fPointRecs[i]->GetZ());
8db76038 510 }
511 fclose(ascii);
512 return;
513}
514//-----------------------------------------------------------------------
515
516void AliITSRiemannFit::ReadPoints(void) {
517 //////////////////////////////////////////////////////////////////////
518 // read the filled array
519 /////////////////////////////////////////////////////////////////////
520 hpsortTrack(fPointRecs,fPoints);
521 for(Int_t i=0;i<fPoints;i++)
522 printf("%d\t%d\t%d\t%f\t%f\t%f\t(%.0f,%.0f,%.0f)\t%.3f\t%s\n",
5ff6c5cc 523 i,fPointRecs[i]->GetLay(),fPointRecs[i]->GetTrack(),
524 fPointRecs[i]->GetX(),fPointRecs[i]->GetY(),
525 fPointRecs[i]->GetZ(),fPointRecs[i]->GetOrigin()->X(),
526 fPointRecs[i]->GetOrigin()->Y(),fPointRecs[i]->GetOrigin()->Z(),
527 fPointRecs[i]->GetPt(),fPointRecs[i]->GetName());
8db76038 528 return;
529}
530//-----------------------------------------------------------------------
531
532Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c,
533 Double_t &x1,Double_t &x2,Double_t &x3){
534 //////////////////////////////////////////////
535 /// Solve cubic equation:
536 /// x^3 + a*x^2 +b*x + c
537 ///
538 /// returns x1 , x2 , x3
539 ////////////////////////////////////////
540
5ff6c5cc 541 Double_t qQ = ((a*a - 3*b)/9);
542 Double_t rR = ((2*a*a*a - 9*a*b +27*c)/54);
8db76038 543 Double_t theta;
5ff6c5cc 544 Double_t aF = -2*sqrt(qQ);
8db76038 545 Double_t g = a/3;
5ff6c5cc 546 Double_t pPI2 = TMath::Pi()*2;
8db76038 547
5ff6c5cc 548 if( rR*rR>qQ*qQ*qQ ) {
8db76038 549 cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<<endl;
550 x1 = 9999999;
551 x2 = 9999999;
552 x3 = 9999999;
553 return 0;
554 }
555
5ff6c5cc 556 theta = TMath::ACos(rR/sqrt(qQ*qQ*qQ));
8db76038 557
5ff6c5cc 558 x1 = (aF*TMath::Cos(theta/3))-g;
559 x2 = (aF*TMath::Cos((theta+pPI2)/3))-g;
560 x3 = (aF*TMath::Cos((theta-pPI2)/3))-g;
8db76038 561
562 return 1;
563}
564//-----------------------------------------------------------------
565
566void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) {
567 ///////////////////////////////////////////////////////////////////////
568 // This function apllies the transformation in the Riemann sphere
569 // for xy plane
570 ///////////////////////////////////////////////////////////////////////
5ff6c5cc 571 Float_t *rR = new Float_t[npoints];
572 Float_t *theta = new Float_t[npoints];
573 Float_t pPI2 = 2*TMath::Pi();
8db76038 574 Float_t x=0,y=0,z=0;
575
576 for(Int_t i=0;i<npoints;i++) {
5ff6c5cc 577 rR[i] = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y());
578 theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X());
579 if(theta[i]<0) theta[i]+=pPI2;
580 x = rR[i]*cos(theta[i])/(1+rR[i]*rR[i]);
581 y = rR[i]*sin(theta[i])/(1+rR[i]*rR[i]);
582 z = rR[i]*rR[i]/(1+rR[i]*rR[i]);
8db76038 583 To[i]->SetXYZ(x,y,z);
584 }
5ff6c5cc 585 delete[] rR;
586 delete[] theta;
8db76038 587 return;
588}
589
590
591//---------------------------------------------------------------------
592
593Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t omega,
594 Double_t &thu0, Double_t &thv0, Double_t &phi, TVector2 &zData, TVector3 &zError,
5ff6c5cc 595 Double_t &corrLin){
8db76038 596 ///////////////////////////////////////////////////////////////////////
597 // Fit the points in the (z,s) plane - helix 3rd equation
598 //
599 ///////////////////////////////////////////////////////////////////////
600 Int_t direction=0;
d65f267e 601 //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
602 //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
603 Double_t * z = new Double_t[npoints];
604 Double_t * x = new Double_t[npoints];
605 Double_t * y = new Double_t[npoints];
606 Double_t * s = new Double_t[npoints];
607 Double_t * ez = new Double_t[npoints];
608 Double_t * ex = new Double_t[npoints];
609 Double_t * ey = new Double_t[npoints];
610 Double_t * es = new Double_t[npoints];
8db76038 611 Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
612
613 // Double_t chi=TMath::Pi()/2.0+phi;
614 Double_t chi=-TMath::Pi()-phi;
615 Double_t angold=0.0, tpang=0.0;
616 for(Int_t k = 0; k<npoints; k++) {
617 x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
618 y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
619 z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
620 if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
621 chisquare=9999.99;
622 cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
d65f267e 623 delete [] z;
624 delete [] x;
625 delete [] y;
626 delete [] s;
627 delete [] ez;
628 delete [] ex;
629 delete [] ey;
630 delete [] es;
8db76038 631 return 12;
632 }
633 Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
634 if( (x[k]-thu0)<0 ) {
635 if (ang1*angold<0) {
636 tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
637 ang1=tpang;
638 }
639 }
640 angold=ang1;
641 if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
642 s[k] = (ang1+chi)/omega;
643 es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
644 }
645 if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
646
647 TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez);
5ff6c5cc 648 fitHist->Fit("pol1","qQ");
8db76038 649 z0 = fitHist->GetFunction("pol1")->GetParameter(0);
650 vpar = fitHist->GetFunction("pol1")->GetParameter(1);
651 ez0 = fitHist->GetFunction("pol1")->GetParError(0);
652 evpar = fitHist->GetFunction("pol1")->GetParError(1);
653 chisquare = fitHist->GetFunction("pol1")->GetChisquare();
654 zData.Set(z0,vpar);
655 zError.SetXYZ(ez0,evpar,chisquare);
656
5ff6c5cc 657 Double_t sigmas=0.;
658 Double_t sigmaz=0.;
659 Double_t avs=0.;
660 Double_t avz=0.;
661 Double_t avsz=0.;
8db76038 662
663 for(Int_t j = 0; j < npoints; j++) {
5ff6c5cc 664 avs += s[j];
665 avz += z[j];
666 avsz += s[j]*z[j];
8db76038 667 }
5ff6c5cc 668 avs /= (Double_t)npoints;
669 avz /= (Double_t)npoints;
670 avsz /= (Double_t)npoints;
8db76038 671
672 for(Int_t l = 0; l < npoints; l++) {
5ff6c5cc 673 sigmas += (s[l]-avs)*(s[l]-avs);
674 sigmaz += (z[l]-avz)*(z[l]-avz);
8db76038 675 }
5ff6c5cc 676 sigmas /=(Double_t)npoints;
677 sigmaz /=(Double_t)npoints;
8db76038 678
5ff6c5cc 679 sigmas = sqrt(sigmas);
680 sigmaz = sqrt(sigmaz);
8db76038 681
5ff6c5cc 682 corrLin = (avsz-avs*avz)/(sigmas*sigmaz);
8db76038 683
d65f267e 684 delete [] z;
685 delete [] x;
686 delete [] y;
687 delete [] s;
688 delete [] ez;
689 delete [] ex;
690 delete [] ey;
691 delete [] es;
692
8db76038 693 return 0;
694}
695
696//-------------------------------------------------------------------
088e0b8d 697Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
8db76038 698 Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,Double_t& omega, Double_t& z0,
699 Double_t& vpar,Double_t& chisql, Double_t& fCorrLin,Double_t& fFit,
d65f267e 700 Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
8db76038 701 ///////////////////////////////////////////////////////////////////////
702 // This function finds the helix paramenters
703 // d0 = impact parameter
704 // rho = radius of circle
705 // phi = atan(y0/x0)
706 // for the xy plane
707 // starting from the momentum and the outcome of
708 // the fit on the Riemann sphere (i.e. u0,v0,rho)
709 //
710 // MIND !!!! Here we assume both angular velocities be 1.0 (yes, one-dot-zero !)
711 //
712 //
713 ///////////////////////////////////////////////////////////////////////
714 //
715 // All this stuff relies on this hypothesis !!!
716 //
717// FILE *pout=fopen("chisql.dat","a");
718 Int_t ierr = 0, ierrl=0;
719 omega = 1.0e-2;
720
721 Int_t bitlay[6]={1,1,1,1,1,1};
722 bitlay[0]*=first; bitlay[1]*=second; bitlay[2]*=third; bitlay[3]*=fourth; bitlay[4]*=fifth; bitlay[5]*=sixth;
723 fd0 = -9999; // No phisycs value
724 u0 = -9999.9999; // parameters of helix - strange value...
725 v0 = -9999.9999; // parameters of helix - strange value...
726 rho = -9999.9999; // parameters of helix -unphysical strange value...
5ff6c5cc 727 Int_t pLayer = 0;
8db76038 728 const Char_t* name = 0;
729 Int_t i=0,k=0;
730 Int_t iMAX = 50;
5ff6c5cc 731 Int_t nN = 0;
8db76038 732 Int_t npl[6]={0,0,0,0,0,0};
5ff6c5cc 733 Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
734 Double_t pt = sqrt(Px*Px+Py*Py);
8db76038 735 TVector3 zError;
736 TVector2 zData;
5ff6c5cc 737 Double_t corrLin;
8db76038 738 TVector3 *ori = new TVector3[iMAX];
739 TVector3 **original = new TVector3*[iMAX];
740 TVector3 *rie = new TVector3[iMAX];
5ff6c5cc 741 TVector3 **riemann = new TVector3*[iMAX];
8db76038 742 TVector3 *err = new TVector3[iMAX];
743 TVector3 **errors = new TVector3*[iMAX];
744 TVector3 *linerr = new TVector3[iMAX];
745 TVector3 **linerrors = new TVector3*[iMAX];
5ff6c5cc 746 //PH Double_t weight[iMAX];
747 Double_t * weight = new Double_t[iMAX];
8db76038 748
749 for(i=0;i<iMAX;i++){
750 original[i] = &(ori[i]);
5ff6c5cc 751 riemann[i] = &(rie[i]);
8db76038 752 errors[i] = &(err[i]);
753 linerrors[i] = &(linerr[i]);
754 }
755 for(k =0;k<iMAX;k++) original[k]->SetXYZ(9999,9999,9999);
5ff6c5cc 756 Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
757 a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
8db76038 758 Double_t xbar = 0;
759 Double_t ybar = 0;
760 Double_t zbar = 0;
761 Double_t a,b,c,d; // cubic parameters
762 Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
763 Double_t value = 0.0; // minimum eigenvalue
764 Double_t x1,x2,x3; // eigenvector component
765 Double_t n1,n2,n3,nr= 0;// unit eigenvector
5ff6c5cc 766 Double_t radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm]
767 Double_t sigmaMS = 0;
768 TVector3 vVec,vVecNor;
8db76038 769
770// Select RecPoints belonging to the track
771 for(k =0;k<fPoints;k++){
5ff6c5cc 772 if(fPointRecs[k]->GetTrack()==tracknumber) {
773 name = fPointRecs[k]->GetName();
774 pt = fPointRecs[k]->GetPt();
775 pLayer = fPointRecs[k]->GetLay();
776 Int_t ilay = pLayer-1;
8db76038 777 if(npl[ilay]!=0) continue;
778 if(bitlay[ilay] == 1) {
5ff6c5cc 779 original[nN]->SetXYZ(0.1*fPointRecs[k]->GetX(),0.1*fPointRecs[k]->GetY(),0.1*fPointRecs[k]->GetZ());
780 errors[nN]->SetXYZ(0.1*fPointRecs[k]->GetdX(),0.1*fPointRecs[k]->GetdY(),0.1*fPointRecs[k]->GetdZ());
781 sigmaMS = (radiusdm[pLayer]-radiusdm[0])*0.000724/pP;// beam pipe contribution
782 for(Int_t j=1;j<pLayer;j++) {
783 sigmaMS += (radiusdm[pLayer]-radiusdm[j])*0.00136/pP;
8db76038 784 }
5ff6c5cc 785 weight[nN] = ( 1 + original[nN]->Perp2() )*( 1+ original[nN]->Perp2() )/
786 ( errors[nN]->Perp2() + sigmaMS*sigmaMS );
787 linerrors[nN]->SetXYZ(errors[nN]->X(),errors[nN]->Y(),sqrt(errors[nN]->Z()*errors[nN]->Z()+sigmaMS*sigmaMS));
788 nN++;
8db76038 789 npl[ilay]++;
790 } // end if on layer
791 } //end if track==tracknumber
792 } //end for k
793 //
794 // 6 points, no more, no less
795 //
796 if(original[5]->X() == 9999 || original[6]->X() != 9999)
797 {
5ff6c5cc 798 delete [] weight;
8db76038 799 return 1; // not enough points
800 }
801
802 //
803 //
804 //
805 // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
806 //
807
5ff6c5cc 808 RiemannTransf(nN,original,riemann);
8db76038 809
5ff6c5cc 810 Double_t sumWeights = 0.0; // sum of weights factor
8db76038 811
5ff6c5cc 812 for(Int_t j=0;j<nN;j++){ // mean values for x[i],y[i],z[i]
813 xbar+=weight[j]*riemann[j]->X();
814 ybar+=weight[j]*riemann[j]->Y();
815 zbar+=weight[j]*riemann[j]->Z();
816 sumWeights+=weight[j];
8db76038 817 }
818
5ff6c5cc 819 xbar /= sumWeights;
820 ybar /= sumWeights;
821 zbar /= sumWeights;
8db76038 822
5ff6c5cc 823 for(Int_t j=0;j<nN;j++) { // Calculate the matrix elements
824 a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
825 a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
826 a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
827 a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
828 a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
829 a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
8db76038 830 }
831
5ff6c5cc 832 a11 /= nN;
833 a12 /= nN;
834 a22 /= nN;
835 a23 /= nN;
836 a13 /= nN;
837 a33 /= nN;
838 a21 = a12;
839 a32 = a23;
840 a31 = a13;
8db76038 841
842// ************** Determinant parameters ********************
843// n.b. simplifications done keeping in mind symmetry of A
844//
845 a = 1;
5ff6c5cc 846 b = (-a11-a33-a22);
847 c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
848 d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
8db76038 849
850// ************** Find the 3 eigenvalues *************************
5ff6c5cc 851 Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
8db76038 852
5ff6c5cc 853 if(checkCubic !=1 ){
8db76038 854 printf("Track %d Has no real solution continuing ...\n",tracknumber);
5ff6c5cc 855 delete [] weight;
8db76038 856 return 2;
857 }
858
859// **************** Find the lowest eigenvalue *****************
860 if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
861 if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
862 if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
863
864 // ************ Eigenvector relative to value **************
865 x3 = 1;
5ff6c5cc 866 x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
867 x1 = (value-a33-a32*x2)/a31;
868 vVec.SetXYZ(x1,x2,x3);
869 vVecNor = vVec.Unit();
870 n1 = vVecNor.X();
871 n2 = vVecNor.Y();
872 n3 = vVecNor.Z();
8db76038 873 nr = -n1*xbar-n2*ybar-n3*zbar;
874
875 u0 = -0.5*n1/(nr+n3);
876 v0 = -0.5*n2/(nr+n3);
877 rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
878
879 fFit = 0.0;
880 fFit += 10.*TMath::Abs(sqrt((original[0]->X()-u0)*(original[0]->X()-u0)+(original[0]->Y()-v0)*(original[0]->Y()-v0))-rho);
881 fFit += 10.*TMath::Abs(sqrt((original[1]->X()-u0)*(original[1]->X()-u0)+(original[1]->Y()-v0)*(original[1]->Y()-v0))-rho);
882 fFit += 10.*TMath::Abs(sqrt((original[2]->X()-u0)*(original[2]->X()-u0)+(original[2]->Y()-v0)*(original[2]->Y()-v0))-rho);
883 fFit += 10.*TMath::Abs(sqrt((original[3]->X()-u0)*(original[3]->X()-u0)+(original[3]->Y()-v0)*(original[3]->Y()-v0))-rho);
884 fFit += 10.*TMath::Abs(sqrt((original[4]->X()-u0)*(original[4]->X()-u0)+(original[4]->Y()-v0)*(original[4]->Y()-v0))-rho);
885 fFit += 10.*TMath::Abs(sqrt((original[5]->X()-u0)*(original[5]->X()-u0)+(original[5]->Y()-v0)*(original[5]->Y()-v0))-rho);
886
887 fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
888 fphi = TMath::ATan2(v0,u0);
889
890//**************************************************************************
891// LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
892// strictly linear (no approximation)
893//**************************************************************************
894
895 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
896 // //
897 // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
898 // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
899 // //
900 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
901
902 rho *= 10.0;
903 u0 *= 10.0;
904 v0 *= 10.0;
5ff6c5cc 905 ierrl=FitLinear(nN,original,linerrors,omega,u0,v0,fphi,zData,zError,corrLin);
8db76038 906 chisql=zError.Z();
907// fprintf(pout,"%f \n",chisql);
908 z0=zData.X();
909 vpar=zData.Y();
5ff6c5cc 910 fCorrLin = corrLin;
8db76038 911 ierr = (ierrl > ierr ? ierrl : ierr);
912// fclose(pout);
5ff6c5cc 913 delete [] weight;
8db76038 914 return ierr;
915}
5ff6c5cc 916Int_t AliITSRiemannFit::FitHelix(Int_t NPoints, TVector3** fPointRecs,TVector3** fPointRecErrors,Float_t& f1, Float_t& f2, Float_t& f3) {
8db76038 917
5ff6c5cc 918 ///////////////////////////////////////////////////////////////////////
919 // This function finds the helix parameters
920 // d0 = impact parameter
921 // rho = radius of circle
922 // phi = atan(y0/x0)
923 // for the xy plane
924 // starting from the momentum and the outcome of
925 // the fit on the Riemann sphere (i.e. u0,v0,rho)
926 //
927 // MIND !!!! Here we assume both angular velocities be 1.0e-2 (yes, 0.01 !)
928 //
929 //
930 // Also linear fit in (z,s) is performed, so it's 3-D !
931 // z0 and vpar are calculated (intercept and z-component of velocity, but
932 // in units... you guess.
933 //
934 //
935 // Values calculated in addition:
936 //
937 // - transverse impact parameter fd0
938 // - sum of residuals in (x,y) plane fFit
939 // - chisquare of linear fit chisql
940 // - correlation coefficient fCorrLin
941 //
942 //
943 //
944 //
945 //
946 ///////////////////////////////////////////////////////////////////////
947 //
948 // All this stuff relies on this hypothesis !!!
949 //
950 Int_t ierr = 0, ierrl=0;
951 const Double_t omega = 1.0e-2;
952
953
954
955
956 Double_t fd0 = -9999; // fake values
957 Double_t u0 = -9999.9999; // for eventual
958 Double_t v0 = -9999.9999; // debugging
959 Double_t rho = -9999.9999; //
960 Double_t fphi, fFit, chisql, z0, vpar, fCorrLin;
961
962 //
963 // This info is no more there... to be re-considered... maybe
964 //
965 // Double_t pP = sqrt(Px*Px+Py*Py+Pz*Pz);
966 // Double_t pt = sqrt(Px*Px+Py*Py);
967
968 TVector3 zError;
969 TVector2 zData;
970 Double_t corrLin;
971 TVector3 *ori = new TVector3[NPoints];
972 TVector3 **original = new TVector3*[NPoints];
973 TVector3 *rie = new TVector3[NPoints];
974 TVector3 **riemann = new TVector3*[NPoints];
975 TVector3 *err = new TVector3[NPoints];
976 TVector3 **errors = new TVector3*[NPoints];
977 TVector3 *linerr = new TVector3[NPoints];
978 TVector3 **linerrors = new TVector3*[NPoints];
979 Double_t * weight = new Double_t[NPoints];
980
981 for(Int_t i=0; i<NPoints; i++){
982
983 original[i] = &(ori[i]);
984 riemann[i] = &(rie[i]);
985 errors[i] = &(err[i]);
986 linerrors[i] = &(linerr[i]);
987
988 original[i]->SetXYZ(9999,9999,9999);
989 }
990
991 //
992 // Riemann fit parameters
993 //
994 Double_t a11,a12,a13,a21,a22,a23,a31,a32,a33;
995 a11=0;a12=0;a13=0;a21=0;a22=0;a23=0;a31=0;a32=0;a33=0;
996 Double_t xbar = 0;
997 Double_t ybar = 0;
998 Double_t zbar = 0;
999 //
1000 Double_t a,b,c,d; // cubic parameters
1001 Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
1002 Double_t value = 0.0; // minimum eigenvalue
1003 Double_t x1,x2,x3; // eigenvector component
1004 Double_t n1,n2,n3,nr= 0; // unit eigenvector
1005 TVector3 vVec,vVecNor;
1006
1007 for (Int_t ip=0; ip<NPoints; ip++) {
1008original[ip]->SetXYZ(0.1*fPointRecs[ip]->X(),0.1*fPointRecs[ip]->Y(),0.1*fPointRecs[ip]->Z());
1009
1010errors[ip]->SetXYZ(0.1*fPointRecErrors[ip]->X(),0.1*fPointRecErrors[ip]->Y(),0.1*fPointRecErrors[ip]->Z());
1011 weight[ip] = (1+original[ip]->Perp2())*(1+original[ip]->Perp2())/(errors[ip]->Perp2());
1012 linerrors[ip]->SetXYZ(errors[ip]->X(),errors[ip]->Y(),errors[ip]->Z());
1013 }
1014
1015
1016 //
1017 //
1018 // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
1019 //
1020
1021 RiemannTransf(NPoints,original,riemann);
1022
1023 Double_t sumWeights = 0.0; // sum of weights factor
1024
1025 for(Int_t j=0;j<NPoints;j++){ // mean values for x[i],y[i],z[i]
1026 xbar+=weight[j]*riemann[j]->X();
1027 ybar+=weight[j]*riemann[j]->Y();
1028 zbar+=weight[j]*riemann[j]->Z();
1029 sumWeights+=weight[j];
1030 }
1031
1032 xbar /= sumWeights;
1033 ybar /= sumWeights;
1034 zbar /= sumWeights;
1035
1036 for(Int_t j=0;j<NPoints;j++) { // Calculate the matrix elements
1037 a11 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->X() - xbar);
1038 a12 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Y() - ybar);
1039 a22 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Y() - ybar);
1040 a23 += weight[j]*(riemann[j]->Y() - ybar)*(riemann[j]->Z() - zbar);
1041 a13 += weight[j]*(riemann[j]->X() - xbar)*(riemann[j]->Z() - zbar);
1042 a33 += weight[j]*(riemann[j]->Z() - zbar)*(riemann[j]->Z() - zbar);
1043 }
1044 //
1045 // this doesn't seem to work...
1046 //
1047 // a11 /= sumWeights;
1048 // a12 /= sumWeights;
1049 // a22 /= sumWeights;
1050 // a23 /= sumWeights;
1051 // a13 /= sumWeights;
1052 // a33 /= sumWeights;
1053
1054 a11 /= NPoints;
1055 a12 /= NPoints;
1056 a22 /= NPoints;
1057 a23 /= NPoints;
1058 a13 /= NPoints;
1059 a33 /= NPoints;
1060 a21 = a12;
1061 a32 = a23;
1062 a31 = a13;
1063
1064// ************** Determinant parameters ********************
1065// n.b. simplifications done keeping in mind symmetry of A
1066//
1067 a = 1;
1068 b = (-a11-a33-a22);
1069 c = (a11*(a22+a33)+a33*a22-a12*a21-a13*a31-a23*a32);
1070 d = (a31*a22*a13+(a12*a21-a11*a22)*a33-2.0*a23*a13*a12+a11*a23*a32);
1071
1072// ************** Find the 3 eigenvalues *************************
1073 Int_t checkCubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
1074
1075 if(checkCubic !=1 ){
1076 printf("No real solution. Check data.\n");
1077 delete [] weight;
1078 return 999;
1079 }
1080
1081// **************** Find the lowest eigenvalue *****************
1082 if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
1083 if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
1084 if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
1085
1086 // ************ Eigenvector relative to value **************
1087 x3 = 1;
1088 x2 = (a33*a21-a23*a31-value*a21)/(a22*a31-a32*a21-value*a31);
1089 x1 = (value-a33-a32*x2)/a31;
1090 vVec.SetXYZ(x1,x2,x3);
1091 vVecNor = vVec.Unit();
1092 n1 = vVecNor.X();
1093 n2 = vVecNor.Y();
1094 n3 = vVecNor.Z();
1095 nr = -n1*xbar-n2*ybar-n3*zbar;
1096
1097 u0 = -0.5*n1/(nr+n3);
1098 v0 = -0.5*n2/(nr+n3);
1099 rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
1100
1101
1102 fFit = 0.0;
1103 for (Int_t i=0; i<NPoints; i++) {
1104 fFit += 10.*TMath::Abs(sqrt((original[i]->X()-u0)*(original[i]->X()-u0)+(original[i]->Y()-v0)*(original[i]->Y()-v0))-rho);
1105 }
1106 fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
1107 fphi = TMath::ATan2(v0,u0);
1108
1109//**************************************************************************
1110// LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
1111// strictly linear (no approximation)
1112//**************************************************************************
1113
1114 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1115 // //
1116 // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
1117 // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
1118 // //
1119 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1120
1121 rho *= 10.0;
1122 u0 *= 10.0;
1123 v0 *= 10.0;
1124
1125 ierrl=LinearFit(NPoints,original,linerrors,omega,u0,v0,fphi,zData,zError,corrLin);
1126 if(ierrl==33) return 0;
1127 chisql=zError.Z();
1128// fprintf(pout,"%f \n",chisql);
1129 z0=zData.X();
1130 vpar=zData.Y();
1131 fCorrLin = corrLin;
1132 ierr = (ierrl > ierr ? ierrl : ierr);
1133// fclose(pout);
1134 delete [] weight;
1135
1136 f1=fphi;
1137 f2=vpar/(omega*TMath::Abs(rho));
1138 f3=1/rho;
1139 delete[] ori;
1140 delete[] rie;
1141 delete[] err;
1142 delete[] linerr;
1143 delete[] original;
1144 delete[] riemann;
1145 delete[] errors;
1146 delete[] linerrors;
1147
1148 return 1;
1149
1150
1151}
1152
1153//____________________________________________________________
1154
1155Int_t AliITSRiemannFit::LinearFit(Int_t npoints, TVector3 **input,
1156 TVector3 **errors, Double_t omega,
1157 Double_t &thu0, Double_t &thv0, Double_t &phi,TVector2 &zData, TVector3 &zError,
1158 Double_t &corrLin){
1159 ///////////////////////////////////////////////////////////////////////
1160 // Fit the points in the (z,s) plane - helix 3rd equation
1161 //
1162 ///////////////////////////////////////////////////////////////////////
1163 //By R.Turrisi
1164
1165 Int_t direction=0;
1166 //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
1167 //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
1168 Double_t * z = new Double_t[npoints];
1169 Double_t * x = new Double_t[npoints];
1170 Double_t * y = new Double_t[npoints];
1171 Double_t * s = new Double_t[npoints];
1172 Double_t * ez = new Double_t[npoints];
1173 Double_t * ex = new Double_t[npoints];
1174 Double_t * ey = new Double_t[npoints];
1175 Double_t * es = new Double_t[npoints];
1176 Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
1177
1178
1179 // Double_t chi=TMath::Pi()/2.0+phi;
1180 Double_t chi=-TMath::Pi()-phi;
1181 Double_t angold=0.0, tpang=0.0;
1182 for(Int_t k = 0; k<npoints; k++) {
1183 x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
1184 y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
1185 z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
1186 if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
1187 chisquare=9999.99;
1188 cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
1189 delete [] z;
1190 delete [] x;
1191 delete [] y;
1192 delete [] s;
1193 delete [] ez;
1194 delete [] ex;
1195 delete [] ey;
1196 delete [] es;
1197 return 12;
1198 }
1199 Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
1200 if( (x[k]-thu0)<0 ) {
1201 if (ang1*angold<0) {
1202 tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
1203 ang1=tpang;
1204 }
8db76038 1205 }
5ff6c5cc 1206 angold=ang1;
1207 if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
1208 s[k] = (ang1+chi)/omega;
1209 es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
1210 }
1211 if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
1212
1213 // if(s[0]>-636 && s[0]<-625) return 33;
1214
1215 TGraph* fitHist = new TGraph(npoints,s,z);
1216 TF1* f1 = new TF1("f1",Fitfunction,-100,100,2);
1217
1218 f1->SetParameter(0,1);
1219 f1->SetParameter(1,1);
1220 f1->SetLineColor(2);
1221 fitHist->Fit(f1,"qQ");
1222
1223 z0 = f1->GetParameter(0);
1224 vpar = f1->GetParameter(1);
1225 ez0 = f1->GetParError(0);
1226 evpar= f1->GetParError(1);
1227 chisquare=f1->GetChisquare();
1228 zData.Set(z0,vpar);
1229 zError.SetXYZ(ez0,evpar,chisquare);
1230
1231 Double_t Sigmas=0.;
1232 Double_t Sigmaz=0.;
1233 Double_t Avs=0.;
1234 Double_t Avz=0.;
1235 Double_t Avsz=0.;
1236
1237 for(Int_t j = 0; j < npoints; j++) {
1238 Avs += s[j];
1239 Avz += z[j];
1240 Avsz += s[j]*z[j];
1241 }
1242 Avs /= (Double_t)npoints;
1243 Avz /= (Double_t)npoints;
1244 Avsz /= (Double_t)npoints;
1245
1246 for(Int_t l = 0; l < npoints; l++) {
1247 Sigmas += (s[l]-Avs)*(s[l]-Avs);
1248 Sigmaz += (z[l]-Avz)*(z[l]-Avz);
1249 }
1250 Sigmas /=(Double_t)npoints;
1251 Sigmaz /=(Double_t)npoints;
1252
1253 Sigmas = sqrt(Sigmas);
1254 Sigmaz = sqrt(Sigmaz);
1255
1256 corrLin = (Avsz-Avs*Avz)/(Sigmas*Sigmaz);
1257
1258
1259
1260 delete [] z;
1261 delete [] x;
1262 delete [] y;
1263 delete [] s;
1264 delete [] ez;
1265 delete [] ex;
1266 delete [] ey;
1267 delete [] es;
1268 delete f1; delete fitHist;
1269 return 0;
8db76038 1270}
5ff6c5cc 1271
1272
1273//_______________________________________________________
1274
1275Double_t AliITSRiemannFit::Fitfunction(Double_t *x, Double_t* par){
1276
1277 return par[0]+(*x)*par[1];
1278
1279}
1280