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