ALIFAST/AliFTrackMaker.cxx

   1 //////////////////////////////////////////////////////////////////////////
   2 //                                                                      //
   3 // AliFast TrackMaker class.                                            //
   4 //                                                                      //
   5 //                                                                      //
   6 //////////////////////////////////////////////////////////////////////////
   7 // ---------------------------------------------------------------------//
   8 //                                                                      //
   9 // origin: "res.f" fortran by Karel Safarik which was used to           //
  10 //         calculate the track resolution for TP.                       //
  11 //         Different detectors and material can be selected.            //
  12 //         The basic routines compute information and error matrices    //
  13 //         used for the calculation of momentum resolution.             //
  14 //         see references: ASK KAREL??                                  //
  15 //                                                                      //
  16 // C++ in AliFast framework: Elzbieta Richter-Was and Yiota Foka        //
  17 //                           following general structure od Makers in   //
  18 //                           ATLFast by R. Brun.                        //
  19 //                                                                      //
  20 // purpose: provide a Maker which by using general basic routines of    //
  21 //          "res.f" computes the necessary elements of covariance matrix//
  22 //          for the calculation of Track Resolution.                    //
  23 //          Those elements are the product of the TrackResolMaker and   //
  24 //          are hold in TrackResol class. They are expected to be used  //
  25 //          together with additional information for the calculation of //
  26 //          the smeared momenta.                                        //
  27 //          Additional information necessary for this calculation       //
  28 //          will be provided via classes or functions specific to the   //
  29 //          specific study and/or detectors.                            //
  30 //          One can select the detector and/or material for a specific  //
  31 //          study.                                                      //
  32 //                                                                      //
  33 // starting point: res.f will be initialy partially contained in        //
  34 //                 AliFTrackResolMaker and in AliFDet                   //
  35 //                 It will be reorganised further in the framework of   //
  36 //                 AliFast according to the needs.                      //
  37 //                 Names of variables are kept as in fortran code.      //
  38 //                                                                      //
  39 //////////////////////////////////////////////////////////////////////////
  40
  41
  42 #ifdef WIN32
  43 // there is a bug in the Microsoft VisualC++ compiler
  44 // this class must be compiled with optimization off on Windows
  45 # pragma optimize( "", off )
  46 #endif
  47
  48 #include <TParticle.h>
  49 #include <TFile.h>
  50 #include <TSystem.h>
  51 #include <TRandom.h>
  52 #include <TROOT.h>
  53 #include <TMath.h>
  54 #include <TH1.h>
  55 #include <TH2.h>
  56 #include <TH3.h>
  57
  58 #include "AliFast.h"
  59 //#include "AliFMCMaker.h"
  60 #include "AliFTrackMaker.h"
  61 #include "AliFTrack.h"
  62 #include "AliFDet.h"
  63 #include "AliMC.h"
  64
  65 const Double_t kPi       = TMath::Pi();
  66 const Double_t k2Pi      = 2*kPi;
  67 const Double_t kPiHalf   = kPi/2.;
  68 extern  AliFast * gAliFast;
  69 ClassImp(AliFTrackMaker)
  70
  71 //_____________________________________________________________________________
  72 AliFTrackMaker::AliFTrackMaker()
  73 {
  74    fNTracks = 0;
  75    fResID1Test = 0;
  76    fResID2Test = 0;
  77    fResID3Test = 0;
  78    fResID4Test = 0;
  79    fResID5Test = 0;
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliFTrackMaker::AliFTrackMaker(const char *name, const char *title)
  84                  :AliFMaker(name,title)
  85 {
  86 //    Default Setters for tracks
  87
  88    fFruits     = new TClonesArray("AliFTrack",100, kFALSE);
  89    fBranchName = "Tracks";
  90    fNTracks    = 0;
  91 // Please, how to do this optionally ??!!!
  92    Save();
  93 }
  94
  95 //_____________________________________________________________________________
  96 AliFTrackMaker::~AliFTrackMaker()
  97 {
  98    //dummy
  99 }
 100
 101 //_____________________________________________________________________________
 102 AliFTrack *AliFTrackMaker::AddTrack(Int_t code, Double_t charge,
 103                                     Double_t pT, Double_t eta,Double_t phi,
 104                                     Double_t v11, Double_t v22, Double_t v33,
 105                                     Double_t v12, Double_t v13, Double_t v23, Int_t iFlag)
 106 {
 107 //            Add a new track to the list of tracks
 108
 109  //Note the use of the "new with placement" to create a new track object.
 110  //This complex "new" works in the following way:
 111  //   tracks[i] is the value of the pointer for track number i in the TClonesArray
 112  //   if it is zero, then a new track must be generated. This typically
 113  //   will happen only at the first events
 114  //   If it is not zero, then the already existing object is overwritten
 115  //   by the new track parameters.
 116  // This technique should save a huge amount of time otherwise spent
 117  // in the operators new and delete.
 118
 119    TClonesArray &tracks = *(TClonesArray*)fFruits;
 120    return new(tracks[fNTracks++]) AliFTrack(code,charge,pT,eta,phi,
 121                                   v11, v22, v33, v12, v13, v23, iFlag);
 122 }
 123
 124 //_____________________________________________________________________________
 125 void AliFTrackMaker::Clear(Option_t *option)
 126 {
 127    //Reset Track Maker
 128
 129    fNTracks = 0;
 130    AliFMaker::Clear(option);
 131 }
 132
 133 //_____________________________________________________________________________
 134 void AliFTrackMaker::Draw(Option_t *)
 135 {
 136 //    Dummy Draw
 137
 138 }
 139
 140 //_____________________________________________________________________________
 141 void AliFTrackMaker::Init()
 142 {
 143   //Create control histograms
 144   if(gAliFast->TestTrack() == 0){
 145
 146      fResID11 = new TH1D("ResID11","Elec: delta(1/pTotal)*pTotal",1000,-0.5,0.5);
 147      fResID12 = new TH1D("ResID12","Elec: delta(lambda)/lambda",1000,-0.01,0.01);
 148      fResID13 = new TH1D("ResID13","Elec: delta(phi)/phi",1000,-0.01,0.01);
 149
 150      fResID21 = new TH1D("ResID21","Pion: delta(1/pTotal)*pTotal",1000,-1.0,1.0);
 151      fResID22 = new TH1D("ResID22","Pion: delta(lambda)/lambda",1000,-1.0,1.0);
 152      fResID23 = new TH1D("ResID23","Pion: delta(phi)/phi",1000,-1.0,1.0);
 153
 154      fResID31 = new TH1D("ResID31","Kaon: delta(1/pTotal)*pTotal",1000,-1.0,1.0);
 155      fResID32 = new TH1D("ResID32","Kaon: delta(lambda)/lambda",1000,-1.0,1.0);
 156      fResID33 = new TH1D("ResID33","Kaon: delta(phi)/phi",1000,-1.0,1.0);
 157
 158      fResID41 = new TH1D("ResID41","Proton: delta(1/pTotal)*pTotal",1000,-1.0,1.0);
 159      fResID42 = new TH1D("ResID42","Proton: delta(lambda)/lambda",1000,-1.0,1.0);
 160      fResID43 = new TH1D("ResID43","Proton: delta(phi)/phi",1000,-1.0,1.0);
 161
 162   }
 163   //Create test histograms for TestJob only
 164   if(gAliFast->TestTrack() == 1){
 165      fResID1Test  = new TH1D("ResID1Test","histogram21 from res.f",1000,0.075,10.075);
 166      fResID2Test  = new TH1D("ResID2Test","histogram21 from res.f",1000,0.075,10.075);
 167      fResID3Test  = new TH1D("ResID3Test","histogram21 from res.f",1000,0.075,10.075);
 168      fResID4Test  = new TH1D("ResID4Test","histogram21 from res.f",1000,0.075,10.075);
 169      fResID5Test  = new TH1D("ResID5Test","histogram21 from res.f",1000,0.075,10.075);
 170   }
 171
 172   //Set particle masses
 173    SetPionMass();
 174    SetKaonMass();
 175    SetElectronMass();
 176    SetProtonMass();
 177
 178   //Switch on/off tracks reconstruction
 179    SetRecTrack();
 180
 181 }
 182
 183 //_____________________________________________________________________________
 184 // Calculate track and its resolution
 185 //_____________________________________________________________________________
 186 void AliFTrackMaker::Make()
 187 {
 188   Double_t v11, v22, v33, v12, v13, v23;
 189   Int_t iFlag;
 190
 191   fNTracks = 0;
 192
 193   // Check if it is a TestJob
 194   if(gAliFast->TestTrack() == 1){
 195      // Run test job
 196      MakeTest(10);
 197   }else{
 198      // Run production job
 199      // Get pointers to Particles arrays and TClonesArray
 200
 201      Int_t idPart, idTrack;
 202      Double_t  charge, pT, eta, phi;
 203      TParticle *part;
 204      Int_t  nparticles = gAlice->GetMCApp()->GetNtrack();
 205      printf("%10s%10d\n","nparticles",nparticles);
 206      for(Int_t ind=0;ind<nparticles;ind++) {
 207        part = gAlice->GetMCApp()->Particle(ind);
 208        idPart  = part->GetPdgCode();
 209        charge  = part->GetPDG()->Charge();
 210        pT      = part->Pt();
 211        eta     = part->Eta();
 212        phi     = part->Phi();
 213        printf("%10s%10d%20.5e%20.5e%20.5e%20.5e\n","Particle",idPart,charge,pT,eta,phi);
 214        // Check convention for tracks reconstruction
 215        idTrack = 0;
 216        if(TMath::Abs(idPart) ==   11)                 idTrack = 1;
 217        if(TMath::Abs(idPart) == 111 || TMath::Abs(idPart) == 211) idTrack = 2;
 218        if(TMath::Abs(idPart) == 311 || TMath::Abs(idPart) == 321) idTrack = 3;
 219        if(TMath::Abs(idPart) == 2212)                 idTrack = 4;
 220
 221        if(idTrack > 0 && fRecTrack > 0){
 222          // Check if track should be reconstructed
 223          if((fRecTrack == 1 && idTrack == 1) ||
 224             (fRecTrack == 2 && idTrack == 2) ||
 225             (fRecTrack == 3 && idTrack == 3) ||
 226             (fRecTrack == 4 && idTrack == 4) ||
 227             fRecTrack == 100 ) {
 228            // Tracks are  reconstructed
 229            ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag);
 230
 231            // Calculate and smear track parameters
 232            Double_t lambda, cosLambda, pTotal,pInverse;
 233            Double_t pInverseSmea, lambdaSmea, phiSmea;
 234            Double_t a1, a2, a3, b2, b3, c3;
 235            Double_t rn1, rn2, rn3;
 236
 237            lambda    = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
 238            cosLambda = TMath::Cos(lambda);
 239            pTotal    = pT/cosLambda;
 240            pInverse  = 1.0/pTotal;
 241
 242            a1  = TMath::Sqrt(v11);
 243            if(a1 == 0.){
 244              a2 = 0;
 245              a3 = 0;
 246            }else{
 247              a2  = v12/a1;
 248              a3  = v13/a1;
 249            }
 250            b2  = TMath::Sqrt(v22-a2*a2);
 251            if(b2 == 0.){
 252              b3 = 0;
 253            }else{
 254              b3  = (v23 - a2*a3)/b2;
 255            }
 256            c3  = TMath::Sqrt(v33 - a3*a3 -b3*b3);
 257            rn1 = gRandom->Gaus(0,1);
 258            rn2 = gRandom->Gaus(0,1);
 259            rn3 = gRandom->Gaus(0,1);
 260
 261            pInverseSmea  = pInverse + a1*rn1;
 262            lambdaSmea    = lambda + a2*rn1 + b2*rn2;
 263            phiSmea       = phi + a3*rn1 + b3*rn2 + c3*rn3;
 264
 265            // Fill control histograms
 266            if(idTrack == 1){
 267              fResID11->Fill((pInverseSmea-pInverse)/pInverse);
 268              fResID12->Fill((lambdaSmea-lambda)/lambda);
 269              fResID13->Fill((phiSmea-phi)/phi);
 270            }
 271            else if(idTrack == 2){
 272              fResID21->Fill((pInverseSmea-pInverse)/pInverse);
 273              fResID22->Fill((lambdaSmea-lambda)/lambda);
 274              fResID23->Fill((phiSmea-phi)/phi);
 275            }
 276            else if(idTrack == 3){
 277              fResID31->Fill((pInverseSmea-pInverse)/pInverse);
 278              fResID32->Fill((lambdaSmea-lambda)/lambda);
 279              fResID33->Fill((phiSmea-phi)/phi);
 280            }
 281            else if(idTrack == 4){
 282              fResID41->Fill((pInverseSmea-pInverse)/pInverse);
 283              fResID42->Fill((lambdaSmea-lambda)/lambda);
 284              fResID43->Fill((phiSmea-phi)/phi);
 285            }
 286          }else{
 287            // Tracks are not reconstructed
 288            v11=0.;
 289            v12=0.;
 290            v13=0.;
 291            v22=0.;
 292            v23=0.;
 293            v33=0.;
 294            iFlag=0;
 295          }
 296          // Store resolution variables  to AliFTrack  ClonesArray
 297          AddTrack(idTrack, charge, pT, eta, phi, v11, v22, v33, v12, v13, v23, iFlag);
 298          printf("%10s%10d%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%10d\n",
 299               "Track",idTrack,charge,pT,eta,phi,v11, v22, v33, v12, v13, v23, iFlag);
 300        }
 301
 302      }
 303   }
 304 }
 305
 306 //_____________________________________________________________________________
 307 void AliFTrackMaker::Finish()
 308 {
 309   // For TestJob only
 310   if(gAliFast->TestTrack() == 1){
 311     /*
 312     // Draw test histograms
 313     TCanvas *c1 = new TCanvas("c1"," ",200,10,600,480);
 314     c1->Divide(2,3);
 315     c1->cd(1);   fResID1Test->Draw();
 316     c1->cd(2);   fResID2Test->Draw();
 317     c1->cd(3);   fResID3Test->Draw();
 318     c1->cd(4);   fResID4Test->Draw();
 319     c1->cd(5);   fResID5Test->Draw();
 320     c1->Update();
 321     // Store TestRes.eps file
 322     c1->Print("TestRes.eps");
 323     */
 324     // Store histograms on file
 325     TFile f2("TestRes.root","RECREATE","Test Res.f");
 326     fResID1Test->Write();
 327     fResID2Test->Write();
 328     fResID3Test->Write();
 329     fResID4Test->Write();
 330     fResID5Test->Write();
 331     f2.Close();
 332   }
 333 }
 334 //_____________________________________________________________________________
 335 void AliFTrackMaker::ErrorMatrix(Int_t idTrack, Double_t pT,  Double_t eta,
 336   Double_t &v11, Double_t &v22, Double_t &v33, Double_t &v12, Double_t &v13, Double_t &v23,
 337   Int_t &iFlag)
 338 {
 339   ///////////////////////////////////////////////
 340   //idTrack      track type            input   //
 341   //pT           transverse mom        input   //
 342   //lambda       deep angle            input   //
 343   //v11,v22,v23  error matrix          output  //
 344   //v12,v13,v23                        output  //
 345   //iFlag                              output  //
 346   ///////////////////////////////////////////////
 347
 348   AliFDet *detector = gAliFast->Detector();
 349   Int_t nDet = detector->NDet();
 350   Int_t nDetActive = detector->NDetActive();
 351   Int_t nTwice = nDetActive + nDetActive;
 352
 353   Double_t rTrack, rTrackInverse, pTotal, pInverse, diffPInverse;
 354   Double_t safety;
 355   Double_t cosLambda, tanLambda, diffLambda;
 356   Double_t rDet;
 357
 358   Double_t hh0[kNMaxDet2][kNMaxDet2], hhi0[kNMaxDet2][kNMaxDet2];
 359   Double_t hh1[kNMaxDet2][kNMaxDet2], hhi1[kNMaxDet2][kNMaxDet2];
 360   Double_t dhhiOverPInverse[kNMaxDet2][kNMaxDet2];
 361   Double_t dhhiOverLambda[kNMaxDet2][kNMaxDet2];
 362   Double_t a1[kNMaxDet2][kNMaxDet2], a2[kNMaxDet2][kNMaxDet2];
 363   Double_t a0PInverse[kNMaxDet2];
 364   Double_t a0Lambda[kNMaxDet2];
 365   Double_t a0Phi[kNMaxDet2];
 366
 367   Double_t vF11, vF12, vF13, vF22, vF23, vF33, d1, d2, d3, det;
 368   Int_t   idet, icyl, im, in;
 369   Double_t phiHalf;
 370   Double_t lambda;
 371
 372   lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
 373   rTrack    = detector->ConstMag()*pT;
 374   safety    = 10.0;
 375   if(2.0*rTrack < (detector->RDet(nDet) + safety)){
 376       iFlag     = 0;
 377       v11 = 0;
 378       v22 = 0;
 379       v33 = 0;
 380       v12 = 0;
 381       v13 = 0;
 382       v23 = 0;
 383       return;
 384   }
 385   iFlag        = 1;
 386   cosLambda    = TMath::Cos(lambda);
 387   pTotal       = pT/cosLambda;
 388   pInverse     = 1.0/pTotal;
 389   diffPInverse = pInverse*1.0e-5;
 390   diffLambda   = 1.0e-4;
 391
 392   // Compute likelihood and derivatives
 393
 394   LogLikelyhood(idTrack, pInverse, lambda);
 395   for(icyl=1; icyl<nTwice+1; icyl++){
 396      for(im=1; im<nTwice+1; im++){
 397       hh0[icyl][im]  = HH(icyl,im);
 398       hhi0[icyl][im] = HHI(icyl,im);
 399      }
 400   }
 401   LogLikelyhood(idTrack, pInverse+diffPInverse,lambda);
 402   for(icyl=1; icyl<nTwice+1; icyl++){
 403      for(im=1; im<nTwice+1; im++){
 404       hh1[icyl][im]  = HH(icyl,im);
 405       hhi1[icyl][im] = HHI(icyl,im);
 406      }
 407   }
 408   for(icyl=1; icyl<nTwice+1; icyl++){
 409      for(im=1; im<icyl+1; im++){
 410         dhhiOverPInverse[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffPInverse;
 411      }
 412   }
 413   LogLikelyhood(idTrack, pInverse, lambda+diffLambda);
 414   for(icyl=1; icyl<nTwice+1; icyl++){
 415      for(im=1; im<nTwice+1; im++){
 416       hh1[icyl][im]  = HH(icyl,im);
 417       hhi1[icyl][im] = HHI(icyl,im);
 418      }
 419   }
 420   for(icyl=1; icyl<nTwice+1; icyl++){
 421      for(im=1; im<icyl+1; im++){
 422         dhhiOverLambda[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffLambda;
 423      }
 424   }
 425
 426   // Compute additional derivatives
 427   rTrackInverse = 1.0/rTrack;
 428   tanLambda    = TMath::Tan(lambda);
 429   icyl = 0;
 430   for(idet=1; idet<nDet+1;idet++){
 431      if(detector->IFlagDet(idet) > 0){
 432         icyl = icyl + 1;
 433         rDet = detector->RDet(idet);
 434         phiHalf = TMath::ASin(0.5*rDet*rTrackInverse);
 435         Double_t rHelp   = rDet /
 436                           (2.0 * TMath::Sqrt(1.0-(0.5 *rDet*rTrackInverse)*
 437                                                  (0.5 *rDet*rTrackInverse)));
 438         a0PInverse[icyl] = - rDet* rHelp
 439                            /(detector->ConstMag()*cosLambda);
 440         a0Lambda[icyl]   = - rDet* rHelp
 441                            * tanLambda * rTrackInverse;
 442         a0Phi[icyl]      =   rDet;
 443         a0PInverse[nDetActive+icyl] = 2.0 * tanLambda
 444                            *rTrack*(rHelp-rTrack*phiHalf)
 445                            /(detector->ConstMag()*cosLambda);
 446         a0Lambda[nDetActive+icyl]   = 2.0 * (  rHelp*tanLambda*tanLambda
 447                                              + rTrack*phiHalf);
 448         a0Phi[nDetActive+icyl] = 0.0 ;
 449     }
 450   }
 451
 452   // Compute information matrix
 453
 454     vF11=0.0;
 455     vF12=0.0;
 456     vF13=0.0;
 457     vF22=0.0;
 458     vF23=0.0;
 459     vF33=0.0;
 460     for(icyl=1; icyl<nTwice+1; icyl++){
 461        d1=0.0;
 462        d2=0.0;
 463        d3=0.0;
 464        for(im=1; im < icyl+1; im++){
 465           d1 = d1 + hhi0[icyl][im]*a0PInverse[im];
 466           d2 = d2 + hhi0[icyl][im]*a0Lambda[im];
 467           d3 = d3 + hhi0[icyl][im]*a0Phi[im];
 468        }
 469        vF11 =vF11 + d1*d1;
 470        vF12 =vF12 + d1*d2;
 471        vF13 =vF13 + d1*d3;
 472        vF22 =vF22 + d2*d2;
 473        vF23 =vF23 + d2*d3;
 474        vF33 =vF33 + d3*d3;
 475     }
 476     for(icyl=1; icyl<nTwice+1; icyl++){
 477        for(im=1; im<icyl+1; im++){
 478           a1[icyl][im] = 0;
 479           a2[icyl][im] = 0;
 480           for(in=im; in<icyl+1;in++){
 481              a1[icyl][im]=a1[icyl][im]+dhhiOverPInverse[icyl][in]*hh0[im][in];
 482              a2[icyl][im]=a2[icyl][im]+dhhiOverLambda[icyl][in]*hh0[im][in];
 483           }
 484           vF11=vF11+a1[icyl][im]*a1[icyl][im];
 485           vF12=vF12+a1[icyl][im]*a2[icyl][im];
 486           vF22=vF22+a2[icyl][im]*a2[icyl][im];
 487        }
 488        vF11=vF11+a1[icyl][icyl]*a1[icyl][icyl];
 489        vF12=vF12+a1[icyl][icyl]*a2[icyl][icyl];
 490        vF22=vF22+a2[icyl][icyl]*a2[icyl][icyl];
 491        }
 492
 493   // Invert information matrix
 494
 495     det=( vF11*vF22 - vF12*vF12 ) *vF33 + (vF12*vF23 - vF13*vF22)*vF13
 496                                         + (vF12*vF13 - vF11*vF23)*vF23;
 497
 498     v11 = (vF22*vF33 - vF23*vF23)/det;
 499     v22 = (vF11*vF33 - vF13*vF13)/det;
 500     v33 = (vF11*vF22 - vF12*vF12)/det;
 501     v12 = (vF13*vF23 - vF12*vF33)/det;
 502     v13 = (vF12*vF23 - vF13*vF22)/det;
 503     v23 = (vF12*vF13 - vF11*vF23)/det;
 504
 505     }
 506 //_____________________________________________________________________________//
 507 void AliFTrackMaker::LogLikelyhood(Int_t idTrack, Double_t pInverse,Double_t lambda)
 508 {
 509   ///////////////////////////////////////////////
 510   //hh           ??                    output  //
 511   //hhi          ??                    output  //
 512   //idTrack       track type           input   //
 513   //pInverse      inverse  momentum    input   //
 514   //lambda       polar angle of track  input   //
 515   ///////////////////////////////////////////////
 516
 517
 518   AliFDet *detector = gAliFast->Detector();
 519   Int_t nDet = detector->NDet();
 520   Int_t nDetActive = detector->NDetActive();
 521   Int_t nTwice = nDetActive + nDetActive;
 522
 523   Double_t    rDet, rDetSQ;
 524   Int_t      idet, icyl, im, imc;
 525   Double_t    cosLambda, tanLambda, pTotal, pT, rTrack, rTrackSQ;
 526   Double_t    beta, overPBeta, rTrackInv, thickCorr, temp1, temp2;
 527   Double_t    partMassSQ;
 528   Double_t    aShelp[kNMaxDet2], dShelp[kNMaxDet2];
 529   Double_t    projXVXT[kNMaxDet2],projYVXT[kNMaxDet2], projZVXT[kNMaxDet2];
 530   Double_t    proj[kNMaxDet2][kNMaxDet2];
 531   Double_t    erroScatt[kNMaxDet2], variance[kNMaxDet2][kNMaxDet2];
 532   Double_t    erroSQ[kNMaxDet2];
 533   Double_t    hh[kNMaxDet2][kNMaxDet2];
 534   Double_t    hhi[kNMaxDet2][kNMaxDet2];
 535   Double_t    errorVX, errorVY, errorVZ;
 536
 537   cosLambda = TMath::Cos(lambda);
 538   tanLambda = TMath::Tan(lambda);
 539   pTotal    = 1.0/pInverse;
 540   pT        = pTotal * cosLambda;
 541   rTrack    = detector->ConstMag() * pTotal * cosLambda;
 542   rTrackSQ  = rTrack * rTrack;
 543   partMassSQ= ParticleMass(idTrack)*ParticleMass(idTrack);
 544   beta      = pTotal / TMath::Sqrt(partMassSQ+pTotal*pTotal);
 545   overPBeta = 1./(pTotal*beta);
 546   rTrackInv = 1./rTrack;
 547   errorVX   = detector->ErrorVertexX();
 548   errorVY   = detector->ErrorVertexY();
 549   errorVZ   = detector->ErrorVertexZ();
 550
 551
 552   erroScatt[0]=0.0;
 553   erroScatt[1]=0.0;
 554   for(idet=1; idet < nDet; idet++){
 555      thickCorr = detector->ThickDet(idet)/TMath::Sqrt(cosLambda*
 556                  TMath::Sqrt(1.0-0.25*(detector->RDetSQ(idet)/rTrackSQ)));
 557      if(detector->IFlagGas(idet) == 0){
 558          thickCorr = thickCorr * (1.3266 + 0.076 * TMath::Log(thickCorr));}
 559      thickCorr = overPBeta * thickCorr;
 560      erroScatt[idet+1]=thickCorr*thickCorr;
 561   }
 562
 563
 564   icyl = 0;
 565   for(idet=1; idet<nDet+1; idet++){
 566     rDet   = detector->RDet(idet);
 567     rDetSQ = rDet*rDet;
 568     dShelp[idet] = TMath::Sqrt(4.0*rTrackSQ-rDetSQ);
 569     aShelp[idet] = TMath::ASin(rDet/(rTrack+rTrack));
 570     if(detector->IFlagDet(idet) > 0) {
 571        icyl = icyl + 1;
 572        projXVXT[icyl] = rDet * rTrackInv;
 573        projXVXT[nDetActive+icyl] = -tanLambda;
 574        temp1 = (rTrackSQ + rTrackSQ - rDetSQ)/dShelp[idet];
 575        temp2 = rDet/dShelp[idet];
 576        projYVXT[icyl] = temp1*rTrackInv;
 577        projYVXT[nDetActive+icyl] = tanLambda * temp2;
 578        projZVXT[icyl] = 0.0;
 579        projZVXT[nDetActive+icyl] = 1.0;
 580        proj[icyl][1] = 0.0;
 581        proj[nDetActive+icyl][0] = 0.0;
 582        proj[nDetActive+icyl][nDet] = 0.0;
 583        proj[icyl][nDet] = 0.0;
 584        for(im=2; im<idet+1; im++){
 585           proj[icyl][im]= (( rDet
 586                            *(rTrackSQ+rTrackSQ-detector->RDetSQ(im-1))
 587                            - detector->RDet(im-1)*temp1*dShelp[im-1])
 588                          /((rTrackSQ + rTrackSQ)*cosLambda));
 589           proj[nDetActive+icyl][im]= 0.5 * detector->RDet(im-1)
 590                                    * rTrackInv
 591                                    * tanLambda * (detector->RDet(im-1)
 592                                    - dShelp[im-1]*temp2)/cosLambda;
 593           proj[nDetActive+icyl][nDet+im]= (rTrack+rTrack) * (aShelp[idet] - aShelp[im-1])
 594              + ( rDet*dShelp[im-1]-detector->RDet(im-1)*dShelp[idet])
 595              *  dShelp[im-1] * tanLambda * tanLambda
 596                / (dShelp[idet] * (rTrack+rTrack));
 597           proj[icyl][nDet+im]= tanLambda
 598                         *  (rDet*detector->RDet(im-1)*dShelp[im-1]
 599                         / (rTrackSQ+rTrackSQ)
 600                         - (rDetSQ + detector->RDetSQ(im-1)
 601                         -  rDetSQ * detector->RDetSQ(im-1)
 602                         / (rTrackSQ+rTrackSQ))
 603                         / dShelp[idet]);
 604        }
 605        for(im=idet+1; im < nDet+1; im++){
 606           proj[icyl][im] = 0.0;
 607           proj[nDetActive+icyl][im] = 0.0;
 608           proj[nDetActive+icyl][nDet+im] = 0.0;
 609           proj[icyl][nDet+im] = 0.0;
 610        }
 611        if(detector->IFlagDet(idet) == 1){
 612           erroSQ[icyl] = detector->ErrorRPhi(idet);
 613           erroSQ[nDetActive+icyl] = detector->ErrorZ(idet);
 614        }else{
 615           TPCResolution(pT, rDet, lambda);
 616           erroSQ[icyl]            = SigmaRPhiSQ();
 617           erroSQ[nDetActive+icyl] = SigmaZSQ();
 618        }
 619        erroSQ[icyl] = erroSQ[icyl] + detector->ErrorR(idet)*temp2*temp2;
 620        erroSQ[nDetActive+icyl] =  erroSQ[nDetActive+icyl]
 621                                 + detector->ErrorR(idet)*tanLambda*tanLambda;
 622     }
 623   }
 624     for(icyl=1; icyl<nTwice+1; icyl++){
 625       for(im=1; im<icyl+1; im++){
 626           variance[icyl][im]=
 627                projXVXT[icyl]*projXVXT[im]*errorVX
 628               +projYVXT[icyl]*projYVXT[im]*errorVY
 629               +projZVXT[icyl]*projZVXT[im]*errorVZ;
 630           for(imc=1; imc<nDet+1; imc++){
 631              variance[icyl][im]=variance[icyl][im]
 632                     +(proj[icyl][imc]*proj[im][imc]
 633                     + proj[icyl][nDet+imc]*proj[im][nDet+imc])
 634                     * erroScatt[imc];
 635           }
 636       }
 637       variance[icyl][icyl] = variance[icyl][icyl]+erroSQ[icyl];
 638     }
 639
 640     for(icyl=1; icyl<nTwice+1; icyl++){
 641       for(im=icyl; im<nTwice+1; im++){
 642           hh[icyl][im]=variance[im][icyl];
 643           for(imc=1; imc<icyl;imc++){
 644                hh[icyl][im]=hh[icyl][im]-hh[imc][icyl]*hh[imc][im];
 645             }
 646             if(im == icyl){
 647                hh[icyl][im] = TMath::Sqrt(hh[icyl][im]);
 648             }  else {
 649                hh[icyl][im] = hh[icyl][im]/hh[icyl][icyl];
 650             }
 651           }
 652     }
 653
 654     for(icyl=1; icyl<nTwice+1; icyl++){
 655         hhi[icyl][icyl] = 1.0 / hh[icyl][icyl];
 656         for(im=1; im<icyl; im++){
 657            hhi[icyl][im] = 0.0;
 658            for(imc=im; imc<icyl; imc++){
 659               hhi[icyl][im] = hhi[icyl][im]-hh[imc][icyl]*hhi[imc][im];
 660            }
 661            hhi[icyl][im] = hhi[icyl][im]*hhi[icyl][icyl];
 662        }
 663     }
 664
 665     for(icyl=1; icyl<nTwice+1; icyl++){
 666       for(im=1; im<nTwice+1; im++){
 667           SetHH(icyl,im,hh[icyl][im]);
 668           SetHHI(icyl,im,hhi[icyl][im]);
 669       }
 670     }
 671
 672
 673 }
 674 //_____________________________________________________________________________
 675 // translation of routine tpc_resolution of res.f
 676 //_____________________________________________________________________________
 677 void AliFTrackMaker::TPCResolution(Double_t pTransv, Double_t radiPad, Double_t lambda)
 678 {
 679   ///////////////////////////////////////////////
 680   //sigmaRPhiSQ  resolution in r-phi   output  //
 681   //sigmaZSQ     resolution in z       output  //
 682   //pTransv      transverse momentum   input   //
 683   //radiPad      radius of pad row     input   //
 684   //lambda       polar angle of track  input   //
 685   //
 686   //units: cm, GeV/c, radian                   //
 687   //parametrisation of TPC resolution          //
 688   //version of 03.07.1995                      //
 689   //source: Marek Kowalski, Karel Safarik      //
 690   ///////////////////////////////////////////////
 691
 692   Double_t aRCoeff=0.41818e-2;
 693   Double_t bRCoeff=0.17460e-4;
 694   Double_t cRCoeff=0.30993e-8;
 695   Double_t dRCoeff=0.41061e-6;
 696   Double_t aZCoeff=0.39610e-2;
 697   Double_t bZCoeff=0.22443e-4;
 698   Double_t cZCoeff=0.51504e-1;
 699
 700   Double_t sigmaRPhiSQ;
 701   Double_t sigmaZSQ;
 702
 703   sigmaRPhiSQ = aRCoeff - bRCoeff * radiPad * TMath::Tan(lambda)+
 704                (cRCoeff * (radiPad/pTransv) + dRCoeff) * radiPad/pTransv;
 705
 706   sigmaZSQ    = aZCoeff - bZCoeff * radiPad * TMath::Tan(lambda)+
 707                 cZCoeff * TMath::Tan(lambda)*TMath::Tan(lambda);
 708
 709   if(sigmaRPhiSQ < 1.0e-6 ) sigmaRPhiSQ = 1.0e-6;
 710   if(sigmaZSQ    < 1.0e-6 ) sigmaZSQ    = 1.0e-6;
 711
 712   sigmaRPhiSQ =   (TMath::Sqrt(sigmaRPhiSQ) + 0.005)
 713                 * (TMath::Sqrt(sigmaRPhiSQ) + 0.005);
 714   sigmaZSQ    =   (TMath::Sqrt(sigmaZSQ)    + 0.005)
 715                 * (TMath::Sqrt(sigmaZSQ)    + 0.005);
 716
 717   SetSigmaRPhiSQ(sigmaRPhiSQ);
 718   SetSigmaZSQ(sigmaZSQ);
 719
 720
 721 }
 722
 723 //_____________________________________________________________________________
 724 // returns the mass given particle ID
 725 //-----------------------------------------------------------------------------
 726 Double_t AliFTrackMaker::ParticleMass(Int_t idTrack)
 727 {
 728   Double_t mass = 0.0;
 729
 730        if(idTrack == 2){ mass = fPionMass;}
 731   else if(idTrack == 3){ mass = fKaonMass;}
 732   else if(idTrack == 1) {mass = fElectronMass;}
 733   else if(idTrack == 4) {mass = fProtonMass;}
 734
 735   return mass;
 736
 737 }
 738
 739 //_____________________________________________________________________________
 740 // returns the rapidity given particle pT, pz
 741 //-----------------------------------------------------------------------------
 742 Double_t AliFTrackMaker::Rapidity(Double_t pt, Double_t pz)
 743 {
 744 //   Compute rapidity
 745
 746    Double_t etalog = TMath::Log((TMath::Sqrt(pt*pt + pz*pz) + TMath::Abs(pz))/pt);
 747    if (pz < 0 ) return -TMath::Abs(etalog);
 748    else         return  TMath::Abs(etalog);
 749 }
 750
 751 //_____________________________________________________________________________
 752 // returns the phi angle given particle px, py
 753 //-----------------------------------------------------------------------------
 754 Double_t AliFTrackMaker::Angle(Double_t x, Double_t y)
 755 {
 756 //   Compute phi angle of particle
 757 // ... this is a copy of function ULANGL
 758 //  .. sign(a,b) = -abs(a) if b <  0
 759 //               =  abs(a) if b >= 0
 760
 761    Double_t angle = 0;
 762    Double_t r = TMath::Sqrt(x*x + y*y);
 763    if (r < 1e-20) return angle;
 764    if (TMath::Abs(x)/r < 0.8) {
 765       angle = TMath::Sign((Double_t)TMath::Abs(TMath::ACos(x/r)), y);
 766    } else {
 767       angle = TMath::ASin(y/r);
 768       if (x < 0 ) {
 769          if(angle >= 0) angle =  kPi - angle;
 770          else           angle = -kPi - angle;
 771       }
 772    }
 773    return angle;
 774 }
 775 //_____________________________________________________________________________
 776 Int_t AliFTrackMaker::Charge(Int_t kf)
 777 {
 778 //...this is copy of function LUCHGE
 779 //...Purpose: to give three times the charge for a particle/parton.
 780
 781   static Int_t kchg[500] = { -1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,-3,0,
 782         0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,3,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,
 783        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 784        2,-1,2,-1,2,3,0,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,
 785        0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,
 786        0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,
 787        0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,
 788        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 789        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 790        3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,
 791        0,3,0,0,0,0,0,0,0,0,-3,0,0,0,0,0,0,0,0,3,0,-3,0,3,-3,0,0,0,3,6,0,
 792        3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,-3,0,3,6,-3,0,3,-3,0,-3,0,3,6,
 793        0,3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 794        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 795        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 796        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
 797
 798 //    extern integer kfcomp_(integer *);
 799   Int_t ipower;
 800   Int_t ret = 0;
 801   Int_t kfa = TMath::Abs(kf);
 802   Int_t kc  = Compress(kfa);
 803
 804 //...Initial values. Simple case of direct readout.
 805   if (kc == 0) {
 806   } else if (kfa <= 100 || kc <= 80 || kc > 100) {
 807      ret = kchg[kc-1];
 808
 809 // ...Construction from quark content for heavy meson, diquark, baryon.
 810   } else if (kfa/1000 % 10 == 0) {
 811         ipower = kfa/100 % 10;
 812         ret = (kchg[kfa / 100 % 10 - 1] - kchg[kfa/10 % 10 - 1])*Int_t(TMath::Power(-1, ipower));
 813   } else if (kfa / 10 % 10 == 0) {
 814         ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1];
 815   } else {
 816         ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1] + kchg[kfa/10 % 10 - 1];
 817   }
 818
 819 // ...Add on correct sign.
 820   if (kf > 0) return ret;
 821   else        return -ret;
 822 }
 823 //_____________________________________________________________________________
 824 Int_t AliFTrackMaker::Compress(Int_t kf)
 825 {
 826 //...this is copy of function LUCOMP
 827 //...Purpose: to compress the standard KF codes for use in mass and decay
 828 //...arrays; also to check whether a given code actually is defined.
 829 //     from BLOCK LUDATA
 830   static Int_t  kchg[500] = { 1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0,
 831         0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 832         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,1,
 833         1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,
 834         1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,
 835         0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,
 836         1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,
 837         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 838         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,
 839         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0,
 840         0,0,1,0,1,1,0,0,0,0,0,0,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,0,0,0,0,1,1,
 841         1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,
 842         1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 843         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 844         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 845         0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
 846   static Int_t kftab[25] = { 211,111,221,311,321,130,310,213,113,223,313,
 847             323,2112,2212,210,2110,2210,110,220,330,440,30443,30553,0,0 };
 848   static Int_t kctab[25] = { 101,111,112,102,103,221,222,121,131,132,122,
 849             123,332,333,281,282,283,284,285,286,287,231,235,0,0 };
 850
 851   Int_t ret = 0;
 852   Int_t kfla, kflb, kflc, kflr, kfls, kfa, ikf;
 853
 854   kfa = TMath::Abs(kf);
 855 //...Simple cases: direct translation or table.
 856   if (kfa == 0 || kfa >= 100000) {
 857      return ret;
 858   } else if (kfa <= 100) {
 859      ret = kfa;
 860      if (kf < 0 && kchg[kfa - 1] == 0)  ret = 0;
 861      return ret;
 862   } else {
 863      for (ikf = 1; ikf <= 23; ++ikf) {
 864         if (kfa == kftab[ikf-1]) {
 865            ret = kctab[ikf-1];
 866            if (kf < 0 && kchg[ret-1] == 0)  ret = 0;
 867            return ret;
 868         }
 869      }
 870   }
 871 // ...Subdivide KF code into constituent pieces.
 872   kfla = kfa / 1000%10;
 873   kflb = kfa / 100%10;
 874   kflc = kfa / 10%10;
 875   kfls = kfa%10;
 876   kflr = kfa / 10000%10;
 877 // ...Mesons.
 878   if (kfa - kflr*10000 < 1000) {
 879      if (kflb == 0 || kflb == 9 || kflc == 0 || kflc == 9) {
 880      } else if (kflb < kflc) {
 881      } else if (kf < 0 && kflb == kflc) {
 882      } else if (kflb == kflc) {
 883         if (kflr == 0 && kfls == 1) {        ret = kflb + 110;
 884         } else if (kflr == 0 && kfls == 3) { ret = kflb + 130;
 885         } else if (kflr == 1 && kfls == 3) { ret = kflb + 150;
 886         } else if (kflr == 1 && kfls == 1) { ret = kflb + 170;
 887         } else if (kflr == 2 && kfls == 3) { ret = kflb + 190;
 888         } else if (kflr == 0 && kfls == 5) { ret = kflb + 210;
 889         }
 890      } else if (kflb <= 5) {
 891         if (kflr == 0 && kfls == 1) {        ret = (kflb-1)*(kflb-2)/2 + 100 + kflc;
 892         } else if (kflr == 0 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 120 + kflc;
 893         } else if (kflr == 1 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 140 + kflc;
 894         } else if (kflr == 1 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 160 + kflc;
 895         } else if (kflr == 2 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 180 + kflc;
 896         } else if (kflr == 0 && kfls == 5) { ret = (kflb-1)*(kflb-2)/2 + 200 + kflc;
 897         }
 898      } else if (kfls == 1 && kflr <= 1 || kfls == 3 && kflr <= 2 || kfls == 5 && kflr == 0) {
 899         ret = kflb + 80;
 900      }
 901 // ...Diquarks.
 902   } else if ((kflr == 0 || kflr == 1) && kflc == 0) {
 903      if (kfls != 1 && kfls != 3) {
 904      } else if (kfla == 9 || kflb == 0 || kflb == 9) {
 905      } else if (kfla < kflb) {
 906      } else if (kfls == 1 && kfla == kflb) {
 907      } else { ret = 90;
 908      }
 909 // ...Spin 1/2 baryons.
 910   } else if (kflr == 0 && kfls == 2) {
 911      if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
 912      } else if (kfla <= kflc || kfla < kflb) {
 913      } else if (kfla >= 6 || kflb >= 4 || kflc >= 4) {
 914          ret = kfla + 80;
 915      } else if (kflb < kflc) {
 916          ret = (kfla+1)*kfla*(kfla-1)/6 + 300 + kflc*(kflc-1)/2 + kflb;
 917      } else {
 918          ret = (kfla+1)*kfla*(kfla-1)/6 + 330 + kflb*(kflb-1)/2 + kflc;
 919      }
 920 // ...Spin 3/2 baryons.
 921   } else if (kflr == 0 && kfls == 4) {
 922      if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
 923      } else if (kfla < kflb || kflb < kflc) {
 924      } else if (kfla >= 6 || kflb >= 4) {
 925          ret = kfla + 80;
 926      } else {
 927          ret = (kfla+1)*kfla*(kfla-1) / 6 + 360 + kflb*(kflb -1) / 2 + kflc;
 928      }
 929   }
 930     return ret;
 931 }
 932
 933 //_____________________________________________________________________________
 934 // TEST JOB: Calculate tracks resolution
 935 //_____________________________________________________________________________
 936 void AliFTrackMaker::MakeTest(Int_t n)
 937 {
 938   Double_t v11, v22, v33, v12, v13, v23;
 939   Int_t iFlag;
 940   Int_t idTrack;
 941   Double_t  pTStart, pT, eta;
 942
 943   Double_t sumDPop,sumDDip,sumDPhi;
 944   Double_t isum,fm;
 945   Double_t pTotal,partMassSQ,beta,lambda;
 946   Double_t dPop,dLop,dDip,dPhi,rho12,rho13,rho23;
 947   Double_t dPPStrag,dPPTot=0;
 948   //  Double_t resol1[1001][11],resol2[10001][11],resol3[1001][11],
 949   //           resol4[1001][11],resol5[10001][11]
 950   Double_t store1[1001],store2[10001],store3[1001],
 951            store4[1001],store5[10001];
 952
 953
 954   idTrack  = 2;
 955   pTStart = 0.07;
 956   for(Int_t istep=1; istep<n; istep++){
 957       if(istep < 100 && istep >  20) istep = istep -1 +  5;
 958       if(istep < 500 && istep > 100) istep = istep -1 + 25;
 959       if(istep <1000 && istep > 500) istep = istep -1 + 100;
 960       pT = pTStart + 0.01*istep;
 961       eta = - 0.044;
 962       sumDPop = 0;
 963       sumDDip = 0;
 964       sumDPhi = 0;
 965       isum    = 0;
 966       for(Int_t in=1; in<11; in++){
 967          eta    = eta + 0.088;
 968          lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
 969          pTotal = pT / TMath::Cos(lambda);
 970          if(idTrack == 1){
 971            dPPStrag = 0.055 /pT;}
 972          else{
 973            partMassSQ = ParticleMass(idTrack)*ParticleMass(idTrack);
 974            beta       = pTotal/ TMath::Sqrt(pTotal*pTotal + partMassSQ);
 975            dPPStrag   = 0.04/(pT*TMath::Power(beta,2.6666666667));
 976          }
 977          ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag);
 978          if(iFlag == 1){
 979             dLop   = TMath::Sqrt(v11);
 980             dDip   = TMath::Sqrt(v22);
 981             dPhi   = TMath::Sqrt(v33);
 982             rho12  = v12/(dLop*dDip);
 983             rho13  = v13/(dLop*dPhi);
 984             rho23  = v23/(dDip*dPhi);
 985             dPop   = 100. *dLop * pTotal;
 986             dDip   = 1000. * dDip;
 987             dPhi   = 1000. * dPhi;
 988             dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);
 989             //            resol1[istep][in] = dPop;
 990             //            resol2[istep][in] = dDip;
 991             //            resol3[istep][in] = dPhi;
 992             //            resol4[istep][in] = dPPTot;
 993             //            resol5[istep][in] = dPPStrag;
 994             sumDPop = sumDPop + dPop;
 995             sumDDip = sumDDip + dDip;
 996             sumDPhi = sumDPhi + dPhi;
 997             isum    = isum + 1;}
 998          else{
 999             printf("%20s %10.5f %10.5f %20s\n","pT,eta",pT,eta,"cannot smear");
1000          }
1001        }
1002        if(isum > 0){
1003          dPop   = sumDPop/isum;
1004          dDip   = sumDDip/isum;
1005          dPhi   = sumDPhi/isum;
1006          dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);}
1007        else{
1008          dPop   = 0;
1009          dDip   = 0;
1010          dPhi   = 0;
1011        }
1012        store1[istep] = dPop;
1013        store2[istep] = dDip;
1014        store3[istep] = dPhi;
1015        store4[istep] = dPPTot;
1016        store5[istep] = dPPStrag;
1017        if(istep > 20 ){
1018           Int_t im = 5;
1019           if(istep > 100) {im =  25;}
1020           if(istep > 500) {im = 100;}
1021           fm = 1./(1.*im);
1022           for(Int_t ist=1; ist<im; ist++){
1023               //               for(Int_t in=1; in < 11; in++){
1024               //                   resol1[istep-im+ist][in] = resol1[istep-im][in]+
1025               //                         ist*fm*(resol1[istep][in]-resol1[istep-im][in]);
1026               //                   resol2[istep-im+ist][in] = resol2[istep-im][in]+
1027               //                         ist*fm*(resol2[istep][in]-resol2[istep-im][in]);
1028               //                   resol3[istep-im+ist][in] = resol3[istep-im][in]+
1029               //                         ist*fm*(resol3[istep][in]-resol3[istep-im][in]);
1030               //                   resol4[istep-im+ist][in] = resol4[istep-im][in]+
1031               //                         ist*fm*(resol4[istep][in]-resol4[istep-im][in]);
1032               //                   resol5[istep-im+ist][in] = resol5[istep-im][in]+
1033               //                         ist*fm*(resol5[istep][in]-resol5[istep-im][in]);
1034               //             }
1035             store1[istep-im+ist]=store1[istep-im]+
1036                                  ist*fm*(store1[istep]-store1[istep-im]);
1037             store2[istep-im+ist]=store2[istep-im]+
1038                                  ist*fm*(store2[istep]-store2[istep-im]);
1039             store3[istep-im+ist]=store3[istep-im]+
1040                                  ist*fm*(store3[istep]-store3[istep-im]);
1041             store4[istep-im+ist]=store4[istep-im]+
1042                                  ist*fm*(store4[istep]-store4[istep-im]);
1043             store5[istep-im+ist]=store5[istep-im]+
1044                                  ist*fm*(store5[istep]-store5[istep-im]);
1045             // Fill control histograms
1046             fResID1Test->Fill(pTStart + 0.01*(istep-im+ist),store1[istep-im+ist]);
1047             fResID2Test->Fill(pTStart + 0.01*(istep-im+ist),store2[istep-im+ist]);
1048             fResID3Test->Fill(pTStart + 0.01*(istep-im+ist),store3[istep-im+ist]);
1049             fResID4Test->Fill(pTStart + 0.01*(istep-im+ist),store4[istep-im+ist]);
1050             fResID5Test->Fill(pTStart + 0.01*(istep-im+ist),store5[istep-im+ist]);
1051           }
1052           printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n",
1053                        "TestTrack:",istep,store1[istep],store2[istep],store3[istep],
1054                                     store4[istep],store5[istep]);
1055        } else {
1056           printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n",
1057                        "TestTrack:",istep,store1[istep],store2[istep],store3[istep],
1058                                     store4[istep],store5[istep]);
1059           fResID1Test->Fill(pT,store1[istep]);
1060           fResID2Test->Fill(pT,store2[istep]);
1061           fResID3Test->Fill(pT,store3[istep]);
1062           fResID4Test->Fill(pT,store4[istep]);
1063           fResID5Test->Fill(pT,store5[istep]);
1064        }
1065   }
1066 }
1067 //_____________________________________________________________________________