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