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