]> 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
negative indexes allowed
[u/mrichter/AliRoot.git] / ITS / AliITSRiemannFit.cxx
CommitLineData
8db76038 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 * *
15 * *
16 * ///////////////////////////////////////////////////////////////////// *
17 * *
18 * This class performs a fast fit of helices going through the <=6 *
19 * points of the ITS, with the goal of studying tracking and *
20 * vertexing performances. *
21 * Generated kinematics is used to take into account different weights *
22 * associated to points in different layers (with different multiple *
23 * scattering-originated errors). *
24 * *
25 * Based on the work by A. Strandlie, R. Fruhwirth *
26 * *
27 * First implementation by N. Bustreo, R. Turrisi - July 2000 *
28 * *
29 * Further modifications by A. Dainese, R. Turrisi *
30 * *
31 * Contact: Rosario Turrisi, rosario.turrisi@pd.infn.it *
32 * *
33 * **************************************************************************/
34//
35//
36// Modified November, 7th 2001 by Rosario Turrisi
37// (rosario.turrisi@pd.infn.it)
38//
39// FitHelix returns different values. 0=ok, >0 =problem
40// void FitLinear -> Int_t FitLinear to give feedback of errors to FitHelix
41//
42//
43// Modified July, 30th 2001 by Rosario Turrisi
44// (rosario.turrisi@pd.infn.it)
45//
46// Fit for z now in (z,s) plane.
47// Returns parameters in order to write the helix equation
48// and find the right phase/initial point.
49//
50// "PROPER WEIGHTS": (1+R^2)^2/(\sigma_x^2 + \sigma_y^2 + \sigma_MS^2)
51//
4ae5bbc4 52#include <Riostream.h>
8db76038 53#include "AliITSRiemannFit.h"
54#include "AliRun.h"
55#include "TClonesArray.h"
56#include "stdio.h"
57#include "stdlib.h"
4ae5bbc4 58#include "Riostream.h"
8db76038 59#include "TH2.h"
60#include "TMath.h"
61#include "TF1.h"
62#include "TGraphErrors.h"
63#include "TMinuit.h"
64#include "TCanvas.h"
65#include "TStyle.h"
66#include "TRandom.h"
67#include "TParticle.h"
68#include "TFile.h"
69#include "AliITSRecPoint.h"
70#include "AliITSgeom.h"
71#include "AliITSmodule.h"
72
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}
93//----------------------------------------------------------------------
94
95AliITSRiemannFit::~AliITSRiemannFit() {
96 ///////////////////////////////////////////////////////////
97 // Default destructor.
98 // if arrays exist delete them. Then set everything to zero.
99 ////////////////////////////////////////////////////////////
100 if(fPointRecs!=0){
101 for(Int_t i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
102 delete[] fPointRecs;
103 } // end if fPointRecs!=0
104 fSizeEvent = 0;
105 fPointRecs = 0;
106 fPoints = 0;
107 fPrimaryTracks = 0;
108 //
109 // test erase
110// fspdi = 0;
111// fspdo = 0;
112 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
113 return;
114}
115//----------------------------------------------------------------------
116
117AliITSRiemannFit::AliITSRiemannFit(Int_t size,Int_t ntracks) {
118 ///////////////////////////////////////////////////////////
119 // Constructor.
120 // Set fSizeEvent to size and fPrimaryTracks to ntracks.
121 // Others to zero.
122 ////////////////////////////////////////////////////////////
123
124 fSizeEvent = size;
125 fPoints = 0;
126 fPrimaryTracks = ntracks;
127 //
128 // test erase
129// fspdi = 0;
130// fspdo = 0;
131 Point_tl *first = new Point_tl[fSizeEvent];
132 Point_tl **PointRecs = new Point_tl*[fSizeEvent];
133 for(Int_t i=0;i<6;i++)fPLay[i] = 0;
134 for(Int_t j=0;j<fSizeEvent;j++) // create an array of struct
135 PointRecs[j] = &(first[j]);
136}
137// ---------------------------------------------------------------------
138
139void FillPoints(Point_tl **Points,Int_t &index,Float_t *xpoint,
140 Float_t *error,
141 TLorentzVector PE,TLorentzVector OT,Int_t *id,
142 Int_t track,const Char_t *name,Int_t code,
143 Float_t phiorigin){
144 ///////////////////////////////////////////////////////////////////////
145 // Fill the structure Point_tl with the proper data
146 //
147 //////////////////////////////////////////////////////////////////////
148 Float_t PI2 = 2.0*TMath::Pi();
149 Float_t phi,r,x,y,z;
150 Int_t i;
151 i = index;
152 x = xpoint[0];
153 y = xpoint[1];
154 z = xpoint[2];
155 r = sqrt(x*x+y*y);
156 phi = TMath::ATan2(y,x);
157 if(phi<0.0) phi += PI2;
158 Points[i]->phi = phi;
159 Points[i]->eta = -0.5*tan(0.5*TMath::ATan2(r,z));
160 Points[i]->fx = x;
161 Points[i]->fy = y;
162 Points[i]->fz = z;
163 Points[i]->fdx = error[0];
164 Points[i]->fdy = error[1];
165 Points[i]->fdz = error[2];
166 Points[i]->fr = r;
167 Points[i]->track = track;
168 Points[i]->lay = id[0],
169 Points[i]->lad = id[1];
170 Points[i]->det = id[2];
171 Points[i]->fMomentum = PE;
172 Points[i]->fOrigin = OT;
173 Points[i]->fPt = sqrt(PE.X()*PE.X()+PE.Y()*PE.Y());
174 Points[i]->fCode = code;
175 Points[i]->fName = name;
176 Points[i]->vertexPhi = phiorigin;
177 index++;
178 return;
179
180}
181// -----------------------------------------------------------------------
182
088e0b8d 183void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS,
8db76038 184 TTree *TR,Int_t nparticles){
185 //////////////////////////////////////////////////////////////////////
186 // Fill the class member fPointRecs with the reconstructed points
187 // Set All other members to the real values
188 //
189 /////////////////////////////////////////////////////////////////////
190 printf("\n ************* Starting Init Points *************\n");
191 TParticle *part;
192 AliITSgeom *gm = (AliITSgeom*)ITS->GetITSgeom();
193 //get pointer to modules array
194 TObjArray *ITSmodules = ITS->GetModules();
195 Int_t nmodules=ITSmodules->GetEntriesFast();
196 printf("nmodules = %d \n",nmodules);
197 // Get the points from points file
198 AliITSmodule *itsModule;
199 Int_t mod,irec;
200 Stat_t nent;
201 AliITSRecPoint *recp;
202 nent=TR->GetEntries();
203 TClonesArray *ITSrec = ITS->RecPoints();
204
205 Int_t TotRP=0;
206 for (mod=0; mod<nmodules; mod++) {
207 itsModule=(AliITSmodule*)ITSmodules->At(mod);
208 ITS->ResetRecPoints();
209 TR->GetEvent(mod);
210 Int_t nrecp = ITSrec->GetEntries();
211 if(!nrecp) continue;
212 TotRP += nrecp;
213 }
214
215 Int_t iMAX = TotRP;
216 fPrimaryTracks = ntracks;
217 fParticles = nparticles;
218 Point_tl *global = new Point_tl[iMAX];
219 fPointRecs = new Point_tl*[iMAX];
220 //
221 // test erase
222// Point_tl *first = new Point_tl[iMAX];
223// Point_tl *second = new Point_tl[iMAX];
224// fspdi = new Point_tl*[iMAX];
225// fspdo = new Point_tl*[iMAX];
226 for(Int_t j=0;j<iMAX;j++) {
227 fPointRecs[j] = &(global[j]);
228 //
229 // test erase
230// fspdi[j] = &(first[j]);
231// fspdo[j] = &(second[j]);
232 }
233
234 Int_t ieta=0,ieta2=0;
235 Int_t i,id[4],idold[4];
236 Int_t track=0;// // track of hit
237 Float_t xpoint[3],error_plus[3],error_minus[3],global_error[3]; // position and error of the point
238 TLorentzVector OT,PE;
239 Float_t locals[3],locals_error[3],locals_plus[3],locals_minus[3]; // local position and local errors
240 Float_t Phi;
241 Int_t code;
242 const char *name;
243 Int_t layer,ladder,detector;
244 Float_t xcluster,zcluster;
245 Int_t num=0,nspdi=0,nspdo=0,nsddi=0,nsddo=0,nssdi=0,nssdo=0;
246
247 for (mod=0; mod<nmodules; mod++) {
248 itsModule=(AliITSmodule*)ITSmodules->At(mod);
249 ITS->ResetRecPoints();
250 TR->GetEvent(mod);
251 Int_t nrecp = ITSrec->GetEntries();
252 if (!nrecp) continue;
253 itsModule->GetID(layer,ladder,detector);
254
255 for (irec=0;irec<nrecp;irec++) {
256 recp = (AliITSRecPoint*)ITSrec->UncheckedAt(irec);
257 track=recp->fTracks[0];
258 if(track <0 ) continue;
259 xcluster=recp->GetX(); // x on cluster
260 zcluster=recp->GetZ(); // z on cluster
261 part = (TParticle*) gAlice->Particle(track);
262 part->ProductionVertex(OT); // set the vertex
263 part->Momentum(PE); // set the vertex momentum
264 name = part->GetName();
265 code = part->GetPdgCode();
266 Phi = part->Phi();
267 id[0]=layer;
268 id[1]=ladder;
269 id[2]=detector;
270 id[3]=irec;
271 locals[0]=xcluster; // x on cluster
272 locals[1]=0.0; // y on cluster
273 locals[2]=zcluster; // z on cluster
274 locals_error[0]=sqrt(recp->GetSigmaX2());
275 locals_error[1]=0.0;
276 locals_error[2]=sqrt(recp->GetSigmaZ2());
277 locals_plus[0]=xcluster+sqrt(recp->GetSigmaX2()); // x on cluster
278 if(layer==1||layer==2) locals_plus[1]=0.0150/2; // y on cluster
279 else if(layer==3||layer==4) locals_plus[1]=0.0280/2; // y on cluster
280 else if(layer==5||layer==6) locals_plus[1]=0.0300/2; // y on cluster
281 locals_plus[2]=zcluster+sqrt(recp->GetSigmaZ2()); // z on cluster
282 locals_minus[0]=xcluster-sqrt(recp->GetSigmaX2()); // x on cluster
283 locals_minus[1]=0.0; // y on cluster
284 locals_minus[2]=zcluster-sqrt(recp->GetSigmaZ2()); // z on cluster
285
286 gm->LtoG(layer,ladder,detector,locals,xpoint);
287 gm->LtoG(layer,ladder,detector,locals_plus,error_plus);
288 gm->LtoG(layer,ladder,detector,locals_minus,error_minus);
289 global_error[0]=0.5*TMath::Abs(error_plus[0]-error_minus[0]);
290 global_error[1]=0.5*TMath::Abs(error_plus[1]-error_minus[1]);
291 global_error[2]=0.5*TMath::Abs(error_plus[2]-error_minus[2]);
292 if(track<ntracks) {
293 if(TMath::Abs(part->Eta())<=1.0) ieta++;
294 if(TMath::Abs(part->Eta())<=0.5) ieta2++;
295 }
296 if(!(id[0]==idold[0]&&id[1]==idold[1]&&
297 id[2]==idold[2]&&id[3]==idold[3])) {
298 FillPoints(fPointRecs,num,xpoint,global_error,PE,OT,id,track,name,code,Phi);
299 //
300 // test erase
301 switch (idold[0]) {
302 case 1:
303 nspdi++;
304 break;
305 case 2:
306 nspdo++;
307 break;
308 case 3:
309 nsddi++;
310 break;
311 case 4:
312 nsddo++;
313 break;
314 case 5:
315 nssdi++;
316 break;
317 case 6:
318 nssdo++;
319 break;
320 }
321// if(idold[0]==1){
322// FillPoints(fspdi,nspdi,xpoint,global_error,PE,OT,id,track,name,code,Phi);
323// }
324// if(idold[0]==2){
325
326// FillPoints(fspdo,nspdo,xpoint,global_error,PE,OT,id,track,name,code,Phi);
327// }
328// if(idold[0]==3){
329// nsddi++;
330// }
331// if(idold[0]==4){
332// nsddo++;
333// }
334// if(idold[0]==5){
335// nssdi++;
336// }
337// if(idold[0]==6){
338// nssdo++;
339// }
340 for(i=0;i<4;i++) idold[i] = id[i];
341 for(i=0;i<3;i++) xpoint[i] = 0.0;
342 } // end if id != idold
343 } // end for irec
344 }// end for mod
345
346 fPoints = num;
347 fSizeEvent = num;
348 fPLay[0] = nspdi ;
349 fPLay[1] = nspdo ;
350 fPLay[2] = nsddi ;
351 fPLay[3] = nsddo ;
352 fPLay[4] = nssdi ;
353 fPLay[5] = nssdo ;
354 printf("%d primary tracks in eta=+-1\n",ieta);
355 printf("%d primary tracks#2 in eta=+-0.5\n",ieta2);
356 printf("\nInitPoints :\n\nPoints on Layer1 : %d on Layer2 : %d\n",nspdi,nspdo);
357 printf("Points on Layer3 : %d on Layer4 : %d\n",nsddi,nsddo);
358 printf("Points on Layer5 : %d on Layer6 : %d\n",nssdi,nssdo);
359 printf("Points on all Layers: %d\n",num);
360 printf("\n ************* Init Points Finished *************\n");
361 return;
362}
363// ------------------------------------------------------------------------
364///////////////////////////////////////////////////////////
365// Functions for sorting the fPointRecs array
366///////////////////////////////////////////////////////////
367Bool_t SortZ(const Point_tl *s1,const Point_tl *s2){
368 // Z sorting function for qsort.
369 Float_t a;
370
371 a = s1->fz - s2->fz;
372 if(a<0.0) return kTRUE;
373 if(a>0.0) return kFALSE;
374 return kFALSE;
375}
376Bool_t SortTrack(const Point_tl *s1,const Point_tl *s2){
377 // track sorting function for qsort.
378 Float_t a;
379
380 a = s1->track - s2->track;
381 if(a<0.0) return kTRUE;
382 if(a>0.0) return kFALSE;
383 return kFALSE;
384}
385void hpsortTrack(Point_tl **ra,Int_t n){
386 Int_t i,ir,j,l;
387 Point_tl *rra;
388
389 if(n<2) return;
390
391 l = ((n-1) >> 1) +1; // divide 2 + 1
392 ir = n-1;
393 for(;;){
394 if(l>0){
395 rra = ra[--l]; // decrement first
396 }else{
397 rra = ra[ir];
398 ra[ir] = ra[0];
399 if(--ir == 0){ // decrement first
400 ra[0] = rra;
401 break;
402 } // if --ra == 0
403 } // end l>0
404 i = l;
405 j = l+1;
406 while(j<=ir){
407 if( j<ir && SortTrack(ra[j],ra[j+1]) ) j++;
408 if( SortTrack(rra,ra[j]) ){
409 ra[i] = ra[j];
410 i = j;
411 j <<= 1; // time 2.
412 }else{
413 break;
414 } // end if func()
415 } // end while
416 ra[i] = rra;
417 } // end for ever
418}
419void hpsortZ(Point_tl **ra,Int_t n){
420 Int_t i,ir,j,l;
421 Point_tl *rra;
422
423 if(n<2) return;
424
425 l = ((n-1) >> 1) +1; // devide 2 + 1
426 ir = n-1;
427 for(;;){
428 if(l>0){
429 rra = ra[--l]; // decrament first
430 }else{
431 rra = ra[ir];
432 ra[ir] = ra[0];
433 if(--ir == 0){ // decrament first
434 ra[0] = rra;
435 break;
436 } // if --ra == 0
437 } // end l>0
438 i = l;
439 j = l+1;
440 while(j<=ir){
441 if( j<ir && SortZ(ra[j],ra[j+1]) ) j++;
442 if( SortZ(rra,ra[j]) ){
443 ra[i] = ra[j];
444 i = j;
445 j <<= 1; // time 2.
446 }else{
447 break;
448 } // end if func()
449 } // end while
450 ra[i] = rra;
451 } // end for ever
452}
453//-----------------------------------------------------------------------
454////////////////////////////////////////////////////////////////////
455// Sorting functions
456///////////////////////////////////////////////////////////////////
457Int_t Partition(Int_t array[],Int_t left,Int_t right){
458 Int_t val = array[left];
459 Int_t lm = left - 1;
460 Int_t rm = right + 1;
461 for(;;) {
462 do
463 rm--;
464 while
465 (array[rm]>val);
466 do
467 lm++;
468 while
469 (array[lm]<val);
470 if(lm<rm){
471 Int_t tempr = array[rm];
472 array[rm]=array[lm];
473 array[lm]=tempr;
474 }
475 else
476 return rm;
477 }
478
479 return 1;
480}
481
482///////////////////////////////////////////////////////////////////////
483
484void AliITSRiemannFit::WritePoints(void) {
485 /////////////////////////////////////////////////////////////////////
486 // write the data in a file (temporary ascii)
487 /////////////////////////////////////////////////////////////////////
488 FILE *ascii= fopen("AsciiPoints.dat","w");
489 for(Int_t i=0;i<fPoints;i++) {
490 fprintf(ascii,"%d\t%d\t%f\t%f\t%f\n",fPointRecs[i]->lay,
491 fPointRecs[i]->track,fPointRecs[i]->fx,fPointRecs[i]->fy,
492 fPointRecs[i]->fz);
493 }
494 fclose(ascii);
495 return;
496}
497//-----------------------------------------------------------------------
498
499void AliITSRiemannFit::ReadPoints(void) {
500 //////////////////////////////////////////////////////////////////////
501 // read the filled array
502 /////////////////////////////////////////////////////////////////////
503 hpsortTrack(fPointRecs,fPoints);
504 for(Int_t i=0;i<fPoints;i++)
505 printf("%d\t%d\t%d\t%f\t%f\t%f\t(%.0f,%.0f,%.0f)\t%.3f\t%s\n",
506 i,fPointRecs[i]->lay,fPointRecs[i]->track,fPointRecs[i]->fx,
507 fPointRecs[i]->fy,fPointRecs[i]->fz,fPointRecs[i]->fOrigin.X(),
508 fPointRecs[i]->fOrigin.Y(),fPointRecs[i]->fOrigin.Z(),
509 fPointRecs[i]->fPt,fPointRecs[i]->fName);
510 return;
511}
512//-----------------------------------------------------------------------
513
514Int_t AliITSRiemannFit::SolveCubic(Double_t a,Double_t b,Double_t c,
515 Double_t &x1,Double_t &x2,Double_t &x3){
516 //////////////////////////////////////////////
517 /// Solve cubic equation:
518 /// x^3 + a*x^2 +b*x + c
519 ///
520 /// returns x1 , x2 , x3
521 ////////////////////////////////////////
522
523 Double_t Q = ((a*a - 3*b)/9);
524 Double_t R = ((2*a*a*a - 9*a*b +27*c)/54);
525 Double_t theta;
526 Double_t F = -2*sqrt(Q);
527 Double_t g = a/3;
528 Double_t PI2 = TMath::Pi()*2;
529
530 if( R*R>Q*Q*Q ) {
531 cout<<"\nTrack "<<"Determinant :\n\t\t No Real Solutions !!!\n"<<endl;
532 x1 = 9999999;
533 x2 = 9999999;
534 x3 = 9999999;
535 return 0;
536 }
537
538 theta = TMath::ACos( R/sqrt(Q*Q*Q));
539
540 x1 = (F*TMath::Cos(theta/3))-g;
541 x2 = (F*TMath::Cos((theta+PI2)/3))-g;
542 x3 = (F*TMath::Cos((theta-PI2)/3))-g;
543
544 return 1;
545}
546//-----------------------------------------------------------------
547
548void RiemannTransf(Int_t npoints,TVector3 **From,TVector3 **To) {
549 ///////////////////////////////////////////////////////////////////////
550 // This function apllies the transformation in the Riemann sphere
551 // for xy plane
552 ///////////////////////////////////////////////////////////////////////
553 Float_t *R = new Float_t[npoints];
554 Float_t *Theta = new Float_t[npoints];
555 Float_t PI2 = 2*TMath::Pi();
556 Float_t x=0,y=0,z=0;
557
558 for(Int_t i=0;i<npoints;i++) {
559 R[i] = sqrt(From[i]->X()*From[i]->X()+From[i]->Y()*From[i]->Y());
560 Theta[i] = TMath::ATan2(From[i]->Y(),From[i]->X());
561 if(Theta[i]<0) Theta[i]+=PI2;
562 x = R[i]*cos(Theta[i])/(1+R[i]*R[i]);
563 y = R[i]*sin(Theta[i])/(1+R[i]*R[i]);
564 z = R[i]*R[i]/(1+R[i]*R[i]);
565 To[i]->SetXYZ(x,y,z);
566 }
567 delete[] R;
568 delete[] Theta;
569 return;
570}
571
572
573//---------------------------------------------------------------------
574
575Int_t FitLinear(Int_t npoints, TVector3 **input, TVector3 **errors, Double_t omega,
576 Double_t &thu0, Double_t &thv0, Double_t &phi, TVector2 &zData, TVector3 &zError,
577 Double_t &CorrLin){
578 ///////////////////////////////////////////////////////////////////////
579 // Fit the points in the (z,s) plane - helix 3rd equation
580 //
581 ///////////////////////////////////////////////////////////////////////
582 Int_t direction=0;
d65f267e 583 //PH Double_t z[npoints],x[npoints],y[npoints],s[npoints];
584 //PH Double_t ez[npoints],ex[npoints],ey[npoints],es[npoints];
585 Double_t * z = new Double_t[npoints];
586 Double_t * x = new Double_t[npoints];
587 Double_t * y = new Double_t[npoints];
588 Double_t * s = new Double_t[npoints];
589 Double_t * ez = new Double_t[npoints];
590 Double_t * ex = new Double_t[npoints];
591 Double_t * ey = new Double_t[npoints];
592 Double_t * es = new Double_t[npoints];
8db76038 593 Double_t z0=0.0,vpar=0.0,ez0=0.0,evpar=0.0, chisquare;
594
595 // Double_t chi=TMath::Pi()/2.0+phi;
596 Double_t chi=-TMath::Pi()-phi;
597 Double_t angold=0.0, tpang=0.0;
598 for(Int_t k = 0; k<npoints; k++) {
599 x[k] = 10.0*input[k]->X(); ex[k] = 10.0*errors[k]->X();
600 y[k] = 10.0*input[k]->Y(); ey[k] = 10.0*errors[k]->Y();
601 z[k] = 10.0*input[k]->Z(); ez[k] = 10.0*errors[k]->Z();
602 if(TMath::Abs(x[k]-thu0)<1.0e-5) { // should never happen, nor give troubles...
603 chisquare=9999.99;
604 cerr<<"limit for x-x_0 "<<x[k]<<" "<<thu0<<endl;
d65f267e 605 delete [] z;
606 delete [] x;
607 delete [] y;
608 delete [] s;
609 delete [] ez;
610 delete [] ex;
611 delete [] ey;
612 delete [] es;
8db76038 613 return 12;
614 }
615 Double_t ang1=TMath::ATan2((y[k]-thv0),(x[k]-thu0));
616 if( (x[k]-thu0)<0 ) {
617 if (ang1*angold<0) {
618 tpang=ang1-TMath::Sign(TMath::Pi()*2.0,ang1);
619 ang1=tpang;
620 }
621 }
622 angold=ang1;
623 if (k>0) direction+=(z[k]>z[k-1] ? 1 : -1);
624 s[k] = (ang1+chi)/omega;
625 es[k]=TMath::Sqrt(ey[k]*ey[k]+ex[k]*ex[k]/TMath::Power((x[k]-thu0),4))*TMath::Abs(s[k]);
626 }
627 if ( TMath::Abs(direction) != (npoints-1) ) {return 11;}
628
629 TGraphErrors *fitHist = new TGraphErrors(npoints,s,z,es,ez);
630 fitHist->Fit("pol1","Q");
631 z0 = fitHist->GetFunction("pol1")->GetParameter(0);
632 vpar = fitHist->GetFunction("pol1")->GetParameter(1);
633 ez0 = fitHist->GetFunction("pol1")->GetParError(0);
634 evpar = fitHist->GetFunction("pol1")->GetParError(1);
635 chisquare = fitHist->GetFunction("pol1")->GetChisquare();
636 zData.Set(z0,vpar);
637 zError.SetXYZ(ez0,evpar,chisquare);
638
639 Double_t Sigmas=0.;
640 Double_t Sigmaz=0.;
641 Double_t Avs=0.;
642 Double_t Avz=0.;
643 Double_t Avsz=0.;
644
645 for(Int_t j = 0; j < npoints; j++) {
646 Avs += s[j];
647 Avz += z[j];
648 Avsz += s[j]*z[j];
649 }
650 Avs /= (Double_t)npoints;
651 Avz /= (Double_t)npoints;
652 Avsz /= (Double_t)npoints;
653
654 for(Int_t l = 0; l < npoints; l++) {
655 Sigmas += (s[l]-Avs)*(s[l]-Avs);
656 Sigmaz += (z[l]-Avz)*(z[l]-Avz);
657 }
658 Sigmas /=(Double_t)npoints;
659 Sigmaz /=(Double_t)npoints;
660
661 Sigmas = sqrt(Sigmas);
662 Sigmaz = sqrt(Sigmaz);
663
664 CorrLin = (Avsz-Avs*Avz)/(Sigmas*Sigmaz);
665
d65f267e 666 delete [] z;
667 delete [] x;
668 delete [] y;
669 delete [] s;
670 delete [] ez;
671 delete [] ex;
672 delete [] ey;
673 delete [] es;
674
8db76038 675 return 0;
676}
677
678//-------------------------------------------------------------------
088e0b8d 679Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
8db76038 680 Double_t& fphi,Double_t& u0, Double_t& v0, Double_t& rho,Double_t& omega, Double_t& z0,
681 Double_t& vpar,Double_t& chisql, Double_t& fCorrLin,Double_t& fFit,
d65f267e 682 Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
8db76038 683 ///////////////////////////////////////////////////////////////////////
684 // This function finds the helix paramenters
685 // d0 = impact parameter
686 // rho = radius of circle
687 // phi = atan(y0/x0)
688 // for the xy plane
689 // starting from the momentum and the outcome of
690 // the fit on the Riemann sphere (i.e. u0,v0,rho)
691 //
692 // MIND !!!! Here we assume both angular velocities be 1.0 (yes, one-dot-zero !)
693 //
694 //
695 ///////////////////////////////////////////////////////////////////////
696 //
697 // All this stuff relies on this hypothesis !!!
698 //
699// FILE *pout=fopen("chisql.dat","a");
700 Int_t ierr = 0, ierrl=0;
701 omega = 1.0e-2;
702
703 Int_t bitlay[6]={1,1,1,1,1,1};
704 bitlay[0]*=first; bitlay[1]*=second; bitlay[2]*=third; bitlay[3]*=fourth; bitlay[4]*=fifth; bitlay[5]*=sixth;
705 fd0 = -9999; // No phisycs value
706 u0 = -9999.9999; // parameters of helix - strange value...
707 v0 = -9999.9999; // parameters of helix - strange value...
708 rho = -9999.9999; // parameters of helix -unphysical strange value...
709 Int_t Layer = 0;
710 const Char_t* name = 0;
711 Int_t i=0,k=0;
712 Int_t iMAX = 50;
713 Int_t N = 0;
714 Int_t npl[6]={0,0,0,0,0,0};
715 Double_t P = sqrt(Px*Px+Py*Py+Pz*Pz);
716 Double_t Pt = sqrt(Px*Px+Py*Py);
717 TVector3 zError;
718 TVector2 zData;
719 Double_t CorrLin;
720 TVector3 *ori = new TVector3[iMAX];
721 TVector3 **original = new TVector3*[iMAX];
722 TVector3 *rie = new TVector3[iMAX];
723 TVector3 **Riemann = new TVector3*[iMAX];
724 TVector3 *err = new TVector3[iMAX];
725 TVector3 **errors = new TVector3*[iMAX];
726 TVector3 *linerr = new TVector3[iMAX];
727 TVector3 **linerrors = new TVector3*[iMAX];
d65f267e 728 //PH Double_t Weight[iMAX];
729 Double_t * Weight = new Double_t[iMAX];
8db76038 730
731 for(i=0;i<iMAX;i++){
732 original[i] = &(ori[i]);
733 Riemann[i] = &(rie[i]);
734 errors[i] = &(err[i]);
735 linerrors[i] = &(linerr[i]);
736 }
737 for(k =0;k<iMAX;k++) original[k]->SetXYZ(9999,9999,9999);
738 Double_t A11,A12,A13,A21,A22,A23,A31,A32,A33;
739 A11=0;A12=0;A13=0;A21=0;A22=0;A23=0;A31=0;A32=0;A33=0;
740 Double_t xbar = 0;
741 Double_t ybar = 0;
742 Double_t zbar = 0;
743 Double_t a,b,c,d; // cubic parameters
744 Double_t roots[3]= {0.0,0.0,0.0}; // cubic solutions
745 Double_t value = 0.0; // minimum eigenvalue
746 Double_t x1,x2,x3; // eigenvector component
747 Double_t n1,n2,n3,nr= 0;// unit eigenvector
748 Double_t Radiusdm[7] = {0.3,0.4,0.7,1.49,2.38,3.91,4.36}; // beam pipe and layers radii [dm]
749 Double_t sigma_MS = 0;
750 TVector3 Vec,VecNor;
751
752// Select RecPoints belonging to the track
753 for(k =0;k<fPoints;k++){
754 if(fPointRecs[k]->track==tracknumber) {
755 name = fPointRecs[k]->fName;
756 Pt = fPointRecs[k]->fPt;
757 Layer = fPointRecs[k]->lay;
758 Int_t ilay = Layer-1;
759 if(npl[ilay]!=0) continue;
760 if(bitlay[ilay] == 1) {
761 original[N]->SetXYZ(0.1*fPointRecs[k]->fx,0.1*fPointRecs[k]->fy,0.1*fPointRecs[k]->fz);
762 errors[N]->SetXYZ(0.1*fPointRecs[k]->fdx,0.1*fPointRecs[k]->fdy,0.1*fPointRecs[k]->fdz);
763 sigma_MS = (Radiusdm[Layer]-Radiusdm[0])*0.000724/P;// beam pipe contribution
764 for(Int_t j=1;j<Layer;j++) {
765 sigma_MS += (Radiusdm[Layer]-Radiusdm[j])*0.00136/P;
766 }
767 Weight[N] = ( 1 + original[N]->Perp2() )*( 1+ original[N]->Perp2() )/
768 ( errors[N]->Perp2() + sigma_MS*sigma_MS );
769 linerrors[N]->SetXYZ(errors[N]->X(),errors[N]->Y(),sqrt(errors[N]->Z()*errors[N]->Z()+sigma_MS*sigma_MS));
770 N++;
771 npl[ilay]++;
772 } // end if on layer
773 } //end if track==tracknumber
774 } //end for k
775 //
776 // 6 points, no more, no less
777 //
778 if(original[5]->X() == 9999 || original[6]->X() != 9999)
779 {
d65f267e 780 delete [] Weight;
8db76038 781 return 1; // not enough points
782 }
783
784 //
785 //
786 //
787 // FIT ON THE RIEMANN SPHERE FOR (x,y) PLANE
788 //
789
790 RiemannTransf(N,original,Riemann);
791
792 Double_t Sum_Weights = 0.0; // sum of weights factor
793
794 for(Int_t j=0;j<N;j++){ // mean values for x[i],y[i],z[i]
795 xbar+=Weight[j]*Riemann[j]->X();
796 ybar+=Weight[j]*Riemann[j]->Y();
797 zbar+=Weight[j]*Riemann[j]->Z();
798 Sum_Weights+=Weight[j];
799 }
800
801 xbar /= Sum_Weights;
802 ybar /= Sum_Weights;
803 zbar /= Sum_Weights;
804
805 for(Int_t j=0;j<N;j++) { // Calculate the matrix elements
806 A11 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->X() - xbar);
807 A12 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->Y() - ybar);
808 A22 += Weight[j]*(Riemann[j]->Y() - ybar)*(Riemann[j]->Y() - ybar);
809 A23 += Weight[j]*(Riemann[j]->Y() - ybar)*(Riemann[j]->Z() - zbar);
810 A13 += Weight[j]*(Riemann[j]->X() - xbar)*(Riemann[j]->Z() - zbar);
811 A33 += Weight[j]*(Riemann[j]->Z() - zbar)*(Riemann[j]->Z() - zbar);
812 }
813
814 A11 /= N;
815 A12 /= N;
816 A22 /= N;
817 A23 /= N;
818 A13 /= N;
819 A33 /= N;
820 A21 = A12;
821 A32 = A23;
822 A31 = A13;
823
824// ************** Determinant parameters ********************
825// n.b. simplifications done keeping in mind symmetry of A
826//
827 a = 1;
828 b = (-A11-A33-A22);
829 c = (A11*(A22+A33)+A33*A22-A12*A21-A13*A31-A23*A32);
830 d = (A31*A22*A13+(A12*A21-A11*A22)*A33-2.0*A23*A13*A12+A11*A23*A32);
831
832// ************** Find the 3 eigenvalues *************************
833 Int_t Check_Cubic = SolveCubic(b,c,d,roots[0],roots[1],roots[2]);
834
835 if(Check_Cubic !=1 ){
836 printf("Track %d Has no real solution continuing ...\n",tracknumber);
d65f267e 837 delete [] Weight;
8db76038 838 return 2;
839 }
840
841// **************** Find the lowest eigenvalue *****************
842 if(roots[0]<=roots[1] && roots[0]<=roots[2]) value = roots[0];
843 if(roots[1]<=roots[0] && roots[1]<=roots[2]) value = roots[1];
844 if(roots[2]<=roots[0] && roots[2]<=roots[1]) value = roots[2];
845
846 // ************ Eigenvector relative to value **************
847 x3 = 1;
848 x2 = (A33*A21-A23*A31-value*A21)/(A22*A31-A32*A21-value*A31);
849 x1 = (value-A33-A32*x2)/A31;
850 Vec.SetXYZ(x1,x2,x3);
851 VecNor = Vec.Unit();
852 n1 = VecNor.X();
853 n2 = VecNor.Y();
854 n3 = VecNor.Z();
855 nr = -n1*xbar-n2*ybar-n3*zbar;
856
857 u0 = -0.5*n1/(nr+n3);
858 v0 = -0.5*n2/(nr+n3);
859 rho = sqrt((n1*n1 + n2*n2 -4*nr*(nr+n3))/(4*(nr+n3)*(nr+n3)));
860
861 fFit = 0.0;
862 fFit += 10.*TMath::Abs(sqrt((original[0]->X()-u0)*(original[0]->X()-u0)+(original[0]->Y()-v0)*(original[0]->Y()-v0))-rho);
863 fFit += 10.*TMath::Abs(sqrt((original[1]->X()-u0)*(original[1]->X()-u0)+(original[1]->Y()-v0)*(original[1]->Y()-v0))-rho);
864 fFit += 10.*TMath::Abs(sqrt((original[2]->X()-u0)*(original[2]->X()-u0)+(original[2]->Y()-v0)*(original[2]->Y()-v0))-rho);
865 fFit += 10.*TMath::Abs(sqrt((original[3]->X()-u0)*(original[3]->X()-u0)+(original[3]->Y()-v0)*(original[3]->Y()-v0))-rho);
866 fFit += 10.*TMath::Abs(sqrt((original[4]->X()-u0)*(original[4]->X()-u0)+(original[4]->Y()-v0)*(original[4]->Y()-v0))-rho);
867 fFit += 10.*TMath::Abs(sqrt((original[5]->X()-u0)*(original[5]->X()-u0)+(original[5]->Y()-v0)*(original[5]->Y()-v0))-rho);
868
869 fd0 = 100000.*(TMath::Sqrt(u0*u0+v0*v0)-rho); // transverse impact parameter in microns
870 fphi = TMath::ATan2(v0,u0);
871
872//**************************************************************************
873// LINEAR FIT IN (z,s) PLANE: z = zData.X() + zData.Y()*s
874// strictly linear (no approximation)
875//**************************************************************************
876
877 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
878 // //
879 // REMEMBER, HERE STILL LENGHTS IN DM'S FOR ___INPUT___ BUT zDATA PARAMETERS ARE RETURNED IN CM'S //
880 // rho, u0, v0 parameters converted right now to cm's... it's a mess, I'll take care, sometimes... //
881 // //
882 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
883
884 rho *= 10.0;
885 u0 *= 10.0;
886 v0 *= 10.0;
887 ierrl=FitLinear(N,original,linerrors,omega,u0,v0,fphi,zData,zError,CorrLin);
888 chisql=zError.Z();
889// fprintf(pout,"%f \n",chisql);
890 z0=zData.X();
891 vpar=zData.Y();
892 fCorrLin = CorrLin;
893 ierr = (ierrl > ierr ? ierrl : ierr);
894// fclose(pout);
d65f267e 895 delete [] Weight;
8db76038 896 return ierr;
897}
898
899//-----------------------------------------------------------------------------
900void AliITSRiemannFit::Streamer(TBuffer &lRb){
901////////////////////////////////////////////////////////////////////////
902// The default Streamer function "written by ROOT" doesn't write out
903// the arrays referenced by pointers. Therefore, a specific Streamer function
904// has to be written. This function should not be modified but instead added
905// on to so that older versions can still be read.
906////////////////////////////////////////////////////////////////////////
907 // Stream an object of class AliITSRiemannFit.
908 Int_t i,j,n;
909 Int_t ii,jj;
910 n=20;
911
912 if (lRb.IsReading()) {
913 Version_t lRv = lRb.ReadVersion(); if (lRv) { }
914 TObject::Streamer(lRb);
915 lRb >> fSizeEvent;
916 lRb >> fPrimaryTracks;
917 lRb >> fPoints;
918 for(i=0;i<6;i++) lRb >> fPLay[i];
919 if(fPointRecs!=0){
920 for(i=0;i<fSizeEvent;i++) delete[] fPointRecs[i];
921 delete[] fPointRecs;
922 } // end if fPointRecs!=0
923 fPointRecs = new Point_tl*[fSizeEvent];
924 for(i=0;i<fSizeEvent;i++){
925 fPointRecs[i] = new Point_tl[n];
926 for(j=0;j<n;j++){
927 lRb >> fPointRecs[i][j].lay;
928 lRb >> fPointRecs[i][j].lad;
929 lRb >> fPointRecs[i][j].det;
930 lRb >> fPointRecs[i][j].track;
931 lRb >> fPointRecs[i][j].fx;
932 lRb >> fPointRecs[i][j].fy;
933 lRb >> fPointRecs[i][j].fz;
934 lRb >> fPointRecs[i][j].fr;
935 lRb >> fPointRecs[i][j].fdE;
936 lRb >> fPointRecs[i][j].fdx;
937 lRb >> fPointRecs[i][j].fdy;
938 lRb >> fPointRecs[i][j].fdz;
939 lRb >> fPointRecs[i][j].fPt;
940 lRb >> (Char_t*)fPointRecs[i][j].fName;
941 for (ii=0;ii<4;ii++)
942 lRb << fPointRecs[i][j].fOrigin[ii];
943 for (jj=0;jj<4;jj++)
944 lRb << fPointRecs[i][j].fMomentum[jj];
945 lRb >> fPointRecs[i][j].fCode;
946 lRb >> fPointRecs[i][j].phi;
947 lRb >> fPointRecs[i][j].eta;
948 lRb >> fPointRecs[i][j].vertexPhi;
949 } //end for j
950 } //end for i
951// if(fspdi!=0){
952// for(i=0;i<fSizeEvent/6;i++) delete[] fspdi[i];
953// delete[] fspdi;
954// } // end if fspdi!=0
955// fspdi = new Point_tl*[fSizeEvent/6];
956// for(i=0;i<fSizeEvent/6;i++){
957// fspdi[i] = new Point_tl[n];
958// for(j=0;j<n;j++){
959// lRb >> fspdi[i][j].lay;
960// lRb >> fspdi[i][j].lad;
961// lRb >> fspdi[i][j].det;
962// lRb >> fspdi[i][j].track;
963// lRb >> fspdi[i][j].fx;
964// lRb >> fspdi[i][j].fy;
965// lRb >> fspdi[i][j].fz;
966// lRb >> fspdi[i][j].fr;
967// lRb >> fspdi[i][j].fdE;
968// lRb >> fspdi[i][j].fdx;
969// lRb >> fspdi[i][j].fdy;
970// lRb >> fspdi[i][j].fdz;
971// lRb >> fspdi[i][j].fPt;
972// for (ii=0;ii<4;ii++)
973// lRb << fspdi[i][j].fOrigin[ii];
974// for (jj=0;jj<4;jj++)
975// lRb << fspdi[i][j].fMomentum[jj];
976// lRb >> fspdi[i][j].fCode;
977// lRb >> (Char_t*)fspdi[i][j].fName;
978// lRb >> fspdi[i][j].phi;
979// lRb >> fspdi[i][j].eta;
980// lRb >> fspdi[i][j].vertexPhi;
981// } //end for j
982// } //end for i
983// if(fspdo!=0){
984// for(i=0;i<fSizeEvent/6;i++) delete[] fspdo[i];
985// delete[] fspdo;
986// } // end if fspdo!=0
987// fspdo = new Point_tl*[fSizeEvent/6];
988// for(i=0;i<fSizeEvent/6;i++){
989// fspdo[i] = new Point_tl[n];
990// for(j=0;j<n;j++){
991// lRb >> fspdo[i][j].lay;
992// lRb >> fspdo[i][j].lad;
993// lRb >> fspdo[i][j].det;
994// lRb >> fspdo[i][j].track;
995// lRb >> fspdo[i][j].fx;
996// lRb >> fspdo[i][j].fy;
997// lRb >> fspdo[i][j].fz;
998// lRb >> fspdo[i][j].fr;
999// lRb >> fspdo[i][j].fdE;
1000// lRb >> fspdo[i][j].fdx;
1001// lRb >> fspdo[i][j].fdy;
1002// lRb >> fspdo[i][j].fdz;
1003// lRb >> fspdo[i][j].fPt;
1004// for (ii=0;ii<4;ii++)
1005// lRb << fspdo[i][j].fOrigin[ii];
1006// for (jj=0;jj<4;jj++)
1007// lRb << fspdo[i][j].fMomentum[jj];
1008// lRb >> fspdo[i][j].fCode;
1009// lRb >> (Char_t*)fspdo[i][j].fName;
1010// lRb >> fspdo[i][j].phi;
1011// lRb >> fspdo[i][j].eta;
1012// lRb >> fspdo[i][j].vertexPhi;
1013// } //end for j
1014// } //end for i
1015 } else {
1016 lRb.WriteVersion(AliITSRiemannFit::IsA());
1017 TObject::Streamer(lRb);
1018 lRb << fSizeEvent;
1019 lRb << fPrimaryTracks;
1020 lRb << fPoints;
1021 for(i=0;i<6;i++) lRb >> fPLay[i];
1022 for(i=0;i<fSizeEvent;i++) for(j=0;j<n;j++){
1023 lRb << fPointRecs[i][j].lay;
1024 lRb << fPointRecs[i][j].lad;
1025 lRb << fPointRecs[i][j].det;
1026 lRb << fPointRecs[i][j].track;
1027 lRb << fPointRecs[i][j].fx;
1028 lRb << fPointRecs[i][j].fy;
1029 lRb << fPointRecs[i][j].fz;
1030 lRb << fPointRecs[i][j].fr;
1031 lRb << fPointRecs[i][j].fdE;
1032 lRb << fPointRecs[i][j].fdx;
1033 lRb << fPointRecs[i][j].fdy;
1034 lRb << fPointRecs[i][j].fdz;
1035 lRb << fPointRecs[i][j].fPt;
1036 for (ii=0;ii<4;ii++)
1037 lRb << fPointRecs[i][j].fOrigin[ii];
1038 for (jj=0;jj<4;jj++)
1039 lRb << fPointRecs[i][j].fMomentum[jj];
1040 lRb << fPointRecs[i][j].fCode;
1041 lRb << fPointRecs[i][j].fName;
1042 lRb << fPointRecs[i][j].phi;
1043 lRb << fPointRecs[i][j].eta;
1044 lRb << fPointRecs[i][j].vertexPhi;
1045 }
1046// for(i=0;i<fSizeEvent/6;i++) for(j=0;j<n;j++){
1047// lRb << fspdi[i][j].lay;
1048// lRb << fspdi[i][j].lad;
1049// lRb << fspdi[i][j].det;
1050// lRb << fspdi[i][j].track;
1051// lRb << fspdi[i][j].fx;
1052// lRb << fspdi[i][j].fy;
1053// lRb << fspdi[i][j].fz;
1054// lRb << fspdi[i][j].fr;
1055// lRb << fspdi[i][j].fdE;
1056// lRb << fspdi[i][j].fdx;
1057// lRb << fspdi[i][j].fdy;
1058// lRb << fspdi[i][j].fdz;
1059// lRb << fspdi[i][j].fPt;
1060// for (ii=0;ii<4;ii++)
1061// lRb << fspdi[i][j].fOrigin[ii];
1062// for (jj=0;jj<4;jj++)
1063// lRb << fspdi[i][j].fMomentum[jj];
1064// lRb << fspdi[i][j].fCode;
1065// lRb << fspdi[i][j].fName;
1066// lRb << fspdi[i][j].phi;
1067// lRb << fspdi[i][j].eta;
1068// lRb << fspdi[i][j].vertexPhi;
1069// }
1070// for(i=0;i<fSizeEvent/6;i++) for(j=0;j<n;j++){
1071// lRb << fspdo[i][j].lay;
1072// lRb << fspdo[i][j].lad;
1073// lRb << fspdo[i][j].det;
1074// lRb << fspdo[i][j].track;
1075// lRb << fspdo[i][j].fx;
1076// lRb << fspdo[i][j].fy;
1077// lRb << fspdo[i][j].fz;
1078// lRb << fspdo[i][j].fr;
1079// lRb << fspdo[i][j].fdE;
1080// lRb << fspdo[i][j].fdx;
1081// lRb << fspdo[i][j].fdy;
1082// lRb << fspdo[i][j].fdz;
1083// lRb << fspdo[i][j].fPt;
1084// for (ii=0;ii<4;ii++)
1085// lRb << fspdo[i][j].fOrigin[ii];
1086// for (jj=0;jj<4;jj++)
1087// lRb << fspdo[i][j].fMomentum[jj];
1088// lRb << fspdo[i][j].fCode;
1089// lRb << fspdo[i][j].fName;
1090// lRb << fspdo[i][j].phi;
1091// lRb << fspdo[i][j].eta;
1092// lRb << fspdo[i][j].vertexPhi;
1093// }
1094 } // end if reading
1095}