]>
Commit | Line | Data |
---|---|---|
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 | ||
73 | ClassImp(AliITSRiemannFit) | |
74 | ||
75 | ||
76 | AliITSRiemannFit::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 | ||
95 | AliITSRiemannFit::~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 | ||
117 | AliITSRiemannFit::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 | ||
139 | void 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 | ||
183 | void AliITSRiemannFit::InitPoints(Int_t evnt,Int_t ntracks,AliITS *ITS, | |
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 | /////////////////////////////////////////////////////////// | |
367 | Bool_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 | } | |
376 | Bool_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 | } | |
385 | void 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 | } | |
419 | void 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 | /////////////////////////////////////////////////////////////////// | |
457 | Int_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 | ||
484 | void 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 | ||
499 | void 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 | ||
514 | Int_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 | ||
548 | void 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 | ||
575 | Int_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 | //------------------------------------------------------------------- | |
679 | Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Int_t charge,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0, | |
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 | //----------------------------------------------------------------------------- | |
900 | void 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 | } |