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