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
  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;
 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];
 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;
 605       delete [] z;
 606       delete [] x;
 607       delete [] y;
 608       delete [] s;
 609       delete [] ez;
 610       delete [] ex;
 611       delete [] ey;
 612       delete [] es;
 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
 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
 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,
 682                                    Int_t first,Int_t second,Int_t third,Int_t fourth,Int_t fifth,Int_t sixth) {
 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];
 728   //PH  Double_t Weight[iMAX];
 729   Double_t * Weight = new Double_t[iMAX];
 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     {
 780       delete [] Weight;
 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);
 837     delete [] Weight;
 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);
 895   delete [] Weight;
 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 }