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