MC-dependent part of AliRun extracted in AliMC (F.Carminati)
[u/mrichter/AliRoot.git] / ITS / AliITSRiemannFit.cxx
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 <Riostream.h>
53 #include "AliITSRiemannFit.h"
54 #include "AliRun.h"
55 #include "TClonesArray.h"
56 #include "stdio.h"
57 #include "stdlib.h"
58 #include "Riostream.h"
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 #include "AliMC.h"
73
74 ClassImp(AliITSRiemannFit)
75
76
77 AliITSRiemannFit::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
96 AliITSRiemannFit::~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
118 AliITSRiemannFit::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
140 void 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
184 void AliITSRiemannFit::InitPoints(Int_t ntracks,AliITS *ITS,
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
262       part   = (TParticle*) gAlice->GetMCApp()->Particle(track);    
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 ///////////////////////////////////////////////////////////
368 Bool_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 }
377 Bool_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 }
386 void 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 }
420 void 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 ///////////////////////////////////////////////////////////////////
458 Int_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
485 void 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
500 void 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
515 Int_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
549 void 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
576 Int_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;
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];
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; 
606       delete [] z;
607       delete [] x;
608       delete [] y;
609       delete [] s;
610       delete [] ez;
611       delete [] ex;
612       delete [] ey;
613       delete [] es;
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
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   
676   return 0;
677 }
678
679 //-------------------------------------------------------------------
680 Int_t AliITSRiemannFit::FitHelix(Int_t tracknumber,Double_t Px,Double_t Py,Double_t Pz,Double_t& fd0,
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,
683                                    Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
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];
729   //PH  Double_t Weight[iMAX];
730   Double_t * Weight = new Double_t[iMAX];
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     {
781       delete [] Weight;
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);
838     delete [] Weight;
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);
896   delete [] Weight;
897   return ierr;
898 }
899
900 //-----------------------------------------------------------------------------
901 void 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 }