TOF/AliTOFReconstructioner.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 /* $Id$ */
  17
  18 //_________________________________________________________________________
  19 // Manager class for TOF reconstruction.
  20 //
  21 //
  22 //-- Authors: Bologna-ITEP-Salerno Group
  23 //
  24 // Description: Manager class for TOF reconstruction (derived from TTask)
  25 // Summary of the main methods:
  26 // - extraction of the TPC (assumed to be) reconstructed tracks
  27 //   comment: it has to me moved as soon as possible into a separate
  28 //   class AliTOFTrackReader (K. Safarik suggestion)
  29 // - geometrical propagation of the above tracks till TOF detector
  30 // - matching of the tracks with the TOF signals
  31 //
  32 // Remark: the GEANT3.21 geometry is used during the geometrical propagation
  33 // of the tracks in order to know the current volume reached by the track.
  34 //
  35 //////////////////////////////////////////////////////////////////////////////
  36
  37 #include <Riostream.h>
  38 #include <stdlib.h>
  39
  40 #include <TBenchmark.h>
  41 #include <TClonesArray.h>
  42 #include <TF1.h>
  43 #include <TF2.h>
  44 #include <TFile.h>
  45 #include <TFolder.h>
  46 #include <TGeant3.h>
  47 #include <TNtuple.h>
  48 #include <TParticle.h>
  49 #include <TROOT.h>
  50 #include <TSystem.h>
  51 #include <TTask.h>
  52 #include <TTree.h>
  53 #include <TVirtualMC.h>
  54
  55 #include "AliConst.h"
  56 #include "AliTOFGeometry.h"
  57 #include "AliDetector.h"
  58 #include "AliHeader.h"
  59 #include "AliLoader.h"
  60 #include "AliRun.h"
  61 #include "AliRun.h"
  62 #include "AliRunLoader.h"
  63 #include "AliTOF.h"
  64 #include "AliTOFHitMap.h"
  65 #include "AliTOFPad.h"
  66 #include "AliTOFRecHit.h"
  67 #include "AliTOFReconstructioner.h"
  68 #include "AliTOFSDigit.h"
  69 #include "AliTOFTrack.h"
  70 #include "AliTOFhit.h"
  71 #include "AliTOFv1.h"
  72 #include "AliTOFv2.h"
  73 #include "AliTOFv2FHoles.h"
  74 #include "AliTOFv3.h"
  75 #include "AliTOFv4.h"
  76 #include "AliTOFv4T0.h"
  77 #include "AliMC.h"
  78
  79 // #include "../TPC/AliTPC.h"
  80 // AliTPChit class or somewhere
  81 // this line has to be commented till TPC will provide fPx fPy fPz and fL in
  82
  83 ClassImp(AliTOFReconstructioner)
  84
  85 //____________________________________________________________________________
  86   AliTOFReconstructioner::AliTOFReconstructioner():TTask("AliTOFReconstructioner","")
  87 {
  88   // default ctor
  89   fNevents = 0 ;
  90   foutputfile  = 0;
  91   foutputntuple= 0;
  92   fZnoise  = 0;
  93   ftail    = 0;
  94 }
  95
  96 //____________________________________________________________________________
  97   AliTOFReconstructioner::AliTOFReconstructioner(char* headerFile, Option_t* opt, char *RecFile ):TTask("AliTOFReconstructioner","")
  98 {
  99   //
 100   // ctor
 101   //
 102   fNevents = 0 ;     // Number of events to reconstruct, 0 means all evens in current file
 103   foutputfile  = 0;
 104   foutputntuple= 0;
 105   fZnoise  = 0;
 106   ftail    = 0;
 107
 108   Init(opt);
 109
 110   // create output file
 111   if (RecFile){
 112     foutputfile= new TFile(RecFile,"RECREATE","root file for matching");
 113   } else {
 114     char outFileName[100];
 115     strcpy(outFileName,"match");
 116     strcat(outFileName,headerFile);
 117     foutputfile= new TFile(outFileName,"RECREATE","root file for matching");
 118   }
 119
 120   // initialize the ALIROOT geometry
 121   gAlice->Init();
 122   gAlice->Print();
 123
 124   CreateNTuple();
 125
 126   // add Task to //root/Tasks folder
 127   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
 128   roottasks->Add(this) ;
 129 }
 130 //____________________________________________________________________________
 131 void AliTOFReconstructioner::Init(Option_t* opt)
 132 {
 133   // Initialize the AliTOFReconstructioner setting parameters for
 134   // reconstruction.
 135   // Option values: Pb-Pb for Pb-Pb events
 136   //                pp    for pp    events
 137
 138   // set common parameters
 139   fdbg=1;
 140   fNevents    = 1;
 141   fFirstEvent = 1;
 142   fLastEvent  = 1;
 143   fTimeResolution =0.120;
 144   fpadefficiency  =0.99 ;
 145   fEdgeEffect     = 2   ;
 146   fEdgeTails      = 0   ;
 147   fHparameter     = 0.4 ;
 148   fH2parameter    = 0.15;
 149   fKparameter     = 0.5 ;
 150   fK2parameter    = 0.35;
 151   fEffCenter      = fpadefficiency;
 152   fEffBoundary    = 0.65;
 153   fEff2Boundary   = 0.90;
 154   fEff3Boundary   = 0.08;
 155   fResCenter      = 50. ;
 156   fResBoundary    = 70. ;
 157   fResSlope       = 40. ;
 158   fTimeWalkCenter = 0.  ;
 159   fTimeWalkBoundary=0.  ;
 160   fTimeWalkSlope  = 0.  ;
 161   fTimeDelayFlag  = 1   ;
 162   fPulseHeightSlope=2.0 ;
 163   fTimeDelaySlope =0.060;
 164   // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.);
 165   fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter);
 166   fChargeSmearing=0.0   ;
 167   fLogChargeSmearing=0.13;
 168   fTimeSmearing   =0.022;
 169   fAverageTimeFlag=0    ;
 170   fChargeFactorForMatching=1;
 171   fTrackingEfficiency=1.0; // 100% TPC tracking efficiency assumed
 172   fSigmavsp = 1.        ;
 173   fSigmaZ   = 0.        ;
 174   fSigmarphi= 0.        ;
 175   fSigmap   = 0.        ;
 176   fSigmaPhi = 0.        ;
 177   fSigmaTheta=0.        ;
 178   fField    = 0.2       ;
 179   // fRadLenTPC : 0.2 includes TRD / 0.03 TPC only
 180   fRadLenTPC=0.06        ; // last value
 181   fCorrectionTRD=0.     ;
 182   fLastTPCRow=111       ;
 183   fRadiusvtxBound=50.   ; // expressed in [cm]
 184   fStep     = 0.1       ; // expressed in [cm] step during propagation of the
 185                           // track inside TOF volumes
 186   fMatchingStyle=2      ;
 187   /* previous values default
 188   fMaxPixels=70000      ;
 189   fMaxAllTracks=70000   ;
 190   fMaxTracks=15000      ;
 191   */
 192   fMaxPixels=165000      ;
 193   fMaxAllTracks=500000   ;
 194   fMaxTracks=15000      ;
 195
 196   fMaxTOFHits=35000     ;
 197   fPBound      =0.0     ; // bending effect: P_t=0.3*z*B*R , z particle charge
 198   fNoiseSlope=20.       ;
 199   // set parameters as specified in opt
 200   //pp case
 201   if(strstr(opt,"pp")){
 202   fMaxTestTracks=500    ;
 203   fNoise    = 26.       ;
 204   fNoiseMeanTof= 26.4       ; // to check
 205   }
 206   //Pb-Pb case
 207   if(strstr(opt,"Pb-Pb")){
 208   fMaxTestTracks=20     ;
 209   fNoise    = 9400.     ;
 210   fNoiseMeanTof= 26.4       ;
 211   }
 212 }
 213
 214 //____________________________________________________________________________
 215   AliTOFReconstructioner::~AliTOFReconstructioner()
 216 {
 217   //
 218   // dtor
 219   //
 220
 221   if (foutputfile)
 222     {
 223       delete foutputfile;
 224       foutputfile = 0;
 225     }
 226   if (foutputntuple)
 227     {
 228       delete foutputntuple;
 229       foutputntuple = 0;
 230     }
 231
 232   if (fZnoise)
 233     {
 234       delete fZnoise;
 235       fZnoise = 0;
 236     }
 237
 238   if (ftail)
 239     {
 240       delete ftail;
 241       ftail = 0;
 242     }
 243 }
 244
 245 //____________________________________________________________________________
 246 void AliTOFReconstructioner::CreateNTuple()
 247 {
 248   //
 249   // Create a Ntuple where information about reconstructed charged particles
 250   // (both primaries and secondaries) are stored
 251   // Variables: event ipart imam xvtx yvtx zvtx pxvtx pyvtx pzvtx time leng matc text mext
 252   // Meaning:
 253   // event - event number (0, 1, ...)
 254   // ipart - PDG code of particles
 255   // imam  - PDG code for the parent
 256   // =0 for primary particle
 257   // xvtx  - x-coordinate of the vertex (cm)
 258   // yvtx  - y-coordinate of the vertex (cm)
 259   // zvtx  - z-coordinate of the vertex (cm)
 260   // pxvtx - x-coordinate of the momentum in the vertex (GeV)
 261   // pyvtx - y-coordinate of the momentum in the vertex (GeV)
 262   // pzvtx - z-coordinate of the momentum in the vertex (GeV)
 263   // time  - time of flight from TOF for given track (ps) - TOF time for the
 264   //         first TOF hit of the track
 265   // leng  - track length to the TOF pixel (cm), evaluate as a sum of the
 266   // track length from the track vertex to TPC and the average
 267   // length of the extrapolated track from TPC to TOF.
 268   // for the track without TOF hits leng=-abs(leng)
 269   // matc  - index of the (TPC track) - (TOF pixel) matching
 270   // =0 for tracks which are not tracks for matching, i.e.
 271   // there is not hit on the TPC or Rvxt>200 cm
 272   // >0 for tracks with positive matching procedure:
 273   //   =1 or 2 for non-identified tracks:
 274   //     =1, if the corresponding pixel is not fired,
 275   //     =2, if the corresponding pixel is also matched to the
 276   //         other track,
 277   //   =3 or 4 for identified tracks:
 278   //     =3, if identified with true time,
 279   //     =4, if identified with wrong time.
 280   // <0 for tracks with negative mathing procedure:
 281   //   =-1, if track do not reach the pixel plate (curved in the
 282   //        magnetic field),
 283   //   =-2, if track is out of z-size of the TOF,
 284   //   =-3, if track is or into the RICH hole, or into the PHOS hole, or in the space between the plates,
 285   //   =-4, if track is into the dead space of the TOF.
 286   // text  - time of fligth from the matching procedure = time of the
 287   //         pixel corresponding to the track (ps)
 288   //         =0 for the tracks with matc<=1
 289   // mext  - mass of the track from the matching procedure
 290   //           =p*sqrt(900*(text/leng)**2-1), if 900*(text/leng)**2-1>=0
 291   //           =-p*sqrt(abs(900*(text/leng)**2-1)), if 900*(text/leng)**2-1<0
 292
 293   foutputntuple= new TNtuple("Ntuple","matching","event:ipart:imam:xvtx:yvtx:zvtx:pxvtx:pyvtx:pzvtx:time:leng:matc:text:mext",2000000); // buffersize set for 25 Pb-Pb events
 294 }
 295
 296 //__________________________________________________________________
 297 Double_t TimeWithTailR(Double_t* x, Double_t* par)
 298 {
 299   // sigma - par[0], alpha - par[1], part - par[2]
 300   //  at x<part*sigma - gauss
 301   //  at x>part*sigma - TMath::Exp(-x/alpha)
 302   Float_t xx =x[0];
 303   Double_t f;
 304   if(xx<par[0]*par[2]) {
 305     f = TMath::Exp(-xx*xx/(2*par[0]*par[0]));
 306   } else {
 307     f = TMath::Exp(-(xx-par[0]*par[2])/par[1]-0.5*par[2]*par[2]);
 308   }
 309   return f;
 310 }
 311
 312 //____________________________________________________________________________
 313 void AliTOFReconstructioner::Exec(const char* datafile, Option_t *option)
 314 {
 315   //
 316   // Performs reconstruction for TOF detector
 317   //
 318   gBenchmark->Start("TOFReconstruction");
 319
 320
 321   AliRunLoader *rl = AliRunLoader::Open(datafile);
 322   if (rl == 0x0)
 323    {
 324      Error("Exec","Can not open session for file %s",datafile);
 325      return;
 326    }
 327   // Get AliRun object from file or create it if not on file
 328   rl->LoadgAlice();
 329   gAlice = rl->GetAliRun();
 330
 331   AliTOF* TOF = (AliTOF *) gAlice->GetDetector ("TOF");
 332   AliDetector* TPC = gAlice->GetDetector("TPC");
 333
 334   if (!TOF) {
 335     Error("AliTOFReconstructioner","TOF not found");
 336     delete rl;
 337     return;
 338   }
 339   if (!TPC) {
 340     Error("AliTOFReconstructioner","TPC Detector not found");
 341     delete rl;
 342     return;
 343   }
 344   AliLoader* tpcloader = rl->GetLoader("TPCLoader");
 345   if (tpcloader == 0x0)
 346    {
 347     Error("AliTOFReconstructioner","Can not get TPC Loader from Run Loader.");
 348     delete rl;
 349     return;
 350    }
 351
 352   AliLoader* tofloader = rl->GetLoader("TOFLoader");
 353   if (tofloader == 0x0)
 354    {
 355     Error("AliTOFReconstructioner","Can not get TOF Loader from Run Loader.");
 356     delete rl;
 357     return;
 358    }
 359
 360   if (fEdgeTails) ftail = new TF1("tail",TimeWithTailR,-2,2,3);
 361
 362   if (fNevents == 0) fNevents = rl->GetNumberOfEvents();
 363   // You have to set the number of event with the ad hoc setter
 364   // see testrecon.C
 365   if (rl->GetHeader() == 0x0) rl->LoadHeader();
 366
 367   tofloader->LoadHits();
 368   tpcloader->LoadHits();
 369
 370   for (Int_t ievent = 0; ievent < fNevents; ievent++) { // start loop on events
 371     rl->GetEvent(ievent);
 372     Int_t nparticles= rl->GetHeader()->GetNtrack();
 373     if (nparticles <= 0) return;
 374
 375     TClonesArray* tofhits=0;
 376     TClonesArray* tpchits=0;
 377
 378     if (TOF) tofhits = TOF->Hits();
 379     if (TPC) tpchits = TPC->Hits();
 380
 381     TTree *TH = tofloader->TreeH();
 382     if (!TH) return;
 383     Int_t ntracks    = (Int_t) (TH->GetEntries()); // primary tracks
 384     cout << "number of primary tracked tracks in current event " << ntracks << endl; // number of primary tracked tracks
 385     // array declaration and initialization
 386     // TOF arrays
 387     //    Int_t mapPixels[AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates()][AliTOFGeometry::NStripC()][AliTOFGeometry::NpadZ()*AliTOFGeometry::NpadX()];
 388
 389     Int_t *** mapPixels = new Int_t**[AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates()];
 390     for (Int_t i=0; i<AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates(); i++) mapPixels[i] = new Int_t*[AliTOFGeometry::NStripC()];
 391     for (Int_t i=0; i<AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates(); i++) {
 392       for (Int_t j=0; j<AliTOFGeometry::NStripC(); j++) {
 393         mapPixels[i][j]= new Int_t[AliTOFGeometry::NpadZ()*AliTOFGeometry::NpadX()];
 394       }
 395     }
 396
 397
 398     // initializing the previous array
 399     for (Int_t i=0;i<AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates();i++) {
 400       for (Int_t j=0;j<AliTOFGeometry::NStripC();j++) {
 401         for (Int_t l=0;l<AliTOFGeometry::NpadZ()*AliTOFGeometry::NpadX();l++) {
 402           mapPixels[i][j][l]=0;
 403         }
 404       }
 405     }
 406
 407     Float_t * toftime = new Float_t[fMaxAllTracks];
 408     InitArray(toftime, fMaxAllTracks);
 409     AliTOFPad* pixelArray = new AliTOFPad[fMaxPixels];
 410     Int_t* iTOFpixel        = new Int_t[fMaxAllTracks];
 411     InitArray(iTOFpixel   , fMaxAllTracks);
 412     Int_t* kTOFhitFirst     = new Int_t[fMaxAllTracks];
 413     InitArray(kTOFhitFirst, fMaxAllTracks);
 414     AliTOFRecHit* hitArray  = new AliTOFRecHit[fMaxTOFHits];
 415     Int_t isHitOnFiredPad=0; // index used to fill hitArray (array used to store informations
 416                              // about pads that contains an hit)
 417     Int_t ntotFiredPads=0;   // index used to fill array -> total number of fired pads (at least one time)
 418
 419     // TPC arrays
 420     AliTOFTrack* trackArray = new AliTOFTrack[fMaxTracks];
 421     Int_t * iparticle = new Int_t[fMaxAllTracks];
 422     InitArray(iparticle,fMaxAllTracks);
 423     Int_t * iTrackPt  = new Int_t[fMaxTracks];
 424     InitArray(iTrackPt, fMaxTracks);  // array
 425     Float_t * ptTrack = new Float_t[fMaxTracks];
 426     InitArray( ptTrack, fMaxTracks);  // array for selected track pt
 427     Int_t   ntotTPCtracks=0; // total number of selected TPC tracks
 428
 429
 430     // reading TOF hits
 431     if(TOF) ReadTOFHits(ntracks, TH, tofhits, mapPixels, kTOFhitFirst, pixelArray, iTOFpixel, toftime, hitArray,isHitOnFiredPad,ntotFiredPads);
 432     cout << "isHitOnFiredPad " << isHitOnFiredPad << " for event " << ievent << endl;
 433
 434     // start debug for adding noise
 435     // adding noise
 436     Int_t nHitsNoNoise=isHitOnFiredPad;
 437
 438
 439     if(fNoise) AddNoiseFromOuter(option,mapPixels,pixelArray,hitArray,isHitOnFiredPad,ntotFiredPads);
 440     cout << "ntotFiredPads after adding noise  " << ntotFiredPads   << " for event " << ievent << endl;
 441     // set the hitArray distance to nearest hit
 442     SetMinDistance(hitArray,nHitsNoNoise);
 443
 444     // these lines has to be commented till TPC will provide fPx fPy fPz
 445     // and fL in AliTPChit class
 446     // reading TPC hits
 447     /*
 448     if(TPC) ReadTPCHits(ntracks, TH, tpchits, iTrackPt, iparticle, ptTrack, trackArray,ntotTPCtracks);
 449     */
 450
 451     // geometrical matching
 452     if(TOF && TPC) Matching(trackArray,hitArray,mapPixels,pixelArray,kTOFhitFirst,ntotFiredPads,iTrackPt,iTOFpixel,ntotTPCtracks);
 453
 454     // fill ntuple with reconstructed particles from current event
 455     FillNtuple(ntracks,trackArray,hitArray,pixelArray,iTOFpixel,iparticle,toftime,ntotFiredPads,ntotTPCtracks);
 456
 457
 458     // free used memory
 459     delete [] toftime;
 460     delete [] pixelArray;
 461     delete [] iTOFpixel;
 462     delete [] kTOFhitFirst;
 463     delete [] hitArray;
 464     delete [] trackArray;
 465     delete [] iparticle;
 466     delete [] iTrackPt;
 467     delete [] ptTrack;
 468
 469    for (Int_t i=0; i<AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates(); i++) {
 470       for (Int_t j=0; j<AliTOFGeometry::NStripC(); j++) {
 471         delete [] mapPixels[i][j];
 472       }
 473     }
 474     for (Int_t i=0; i<AliTOFGeometry::NSectors()*AliTOFGeometry::NPlates(); i++) delete [] mapPixels[i];
 475
 476     delete [] mapPixels;
 477
 478   }//event loop
 479
 480   // free used memory for ftail
 481   if (ftail)
 482     {
 483       delete ftail;
 484       ftail = 0;
 485     }
 486
 487   // writing ntuple on output file
 488   foutputfile->cd();
 489   //foutputntuple->Write(0,TObject::kOverwrite);
 490   foutputntuple->Write();
 491   foutputfile->Write();
 492   foutputfile->Close();
 493
 494   gBenchmark->Stop("TOFReconstruction");
 495   cout << "AliTOFReconstructioner:" << endl ;
 496   cout << "   took " << gBenchmark->GetCpuTime("TOFReconstruction") << " seconds in order to make the reconstruction for " <<  fNevents << " events " << endl;
 497   cout <<  gBenchmark->GetCpuTime("TOFReconstruction")/fNevents << " seconds per event " << endl ;
 498   cout << endl ;
 499
 500 }
 501
 502 //__________________________________________________________________
 503 void AliTOFReconstructioner::SetRecFile(char * file )
 504 {
 505   //
 506   // Set the file name for reconstruction output
 507   //
 508   if(!fRecFile.IsNull())
 509     cout << "Changing destination file for TOF reconstruction from " <<(char *)fRecFile.Data() << " to " << file << endl ;
 510   fRecFile=file ;
 511 }
 512 //__________________________________________________________________
 513 void AliTOFReconstructioner::Print(Option_t* /*option*/)const
 514 {
 515   //
 516   // Print reconstruction output file name
 517   //
 518   cout << "------------------- "<< GetName() << " -------------" << endl ;
 519   if(fRecFile.IsNull())
 520     cout << " Writing reconstructed particles to file galice.root "<< endl ;
 521   else
 522     cout << "    Writing reconstructed particle to file  " << (char*) fRecFile.Data() << endl ;
 523
 524 }
 525
 526 //__________________________________________________________________
 527 void AliTOFReconstructioner::PrintParameters()const
 528 {
 529   //
 530   // Print parameters used for reconstruction
 531   //
 532   cout << " ------------------- "<< GetName() << " -------------" << endl ;
 533   cout << " Parameters used for TOF reconstruction " << endl ;
 534   //  Printing the parameters
 535
 536   cout << " Number of events:                        " << fNevents << endl;
 537   cout << " Recostruction from event                 "<< fFirstEvent << "  to event "<< fLastEvent << endl;
 538   cout << " TOF geometry parameters                  " << endl;
 539   cout << " Min. radius of the TOF (cm)              "<< AliTOFGeometry::Rmin() << endl;
 540   cout << " Max. radius of the TOF (cm)              "<< AliTOFGeometry::Rmax() << endl;
 541   cout << " Number of TOF geom. levels               "<< AliTOFGeometry::MaxTOFTree()<< endl;
 542   cout << " Number of TOF sectors                    "<< AliTOFGeometry::NSectors() << endl;
 543   cout << " Number of TOF modules                    "<< AliTOFGeometry::NPlates() << endl;
 544   cout << " Max. Number of strips in a module        "<< AliTOFGeometry::NStripC() << endl;
 545   cout << " Number of pads per strip                 "<< AliTOFGeometry::NpadX()*AliTOFGeometry::NpadZ() << endl;
 546   cout << " Number of strips in central module       "<< AliTOFGeometry::NStripA() << endl;
 547   cout << " Number of strips in intermediate modules "<< AliTOFGeometry::NStripB() << endl;
 548   cout << " Number of strips in outer modules        "<< AliTOFGeometry::NStripC() << endl;
 549   cout << " Number of MRPC in x strip direction      "<< AliTOFGeometry::NpadX()<< endl;
 550   cout << " Size of MRPC (cm) along X                "<< AliTOFGeometry::XPad()<< endl;
 551   cout << " Number of MRPC in z strip direction      "<< AliTOFGeometry::NpadZ()<<endl;
 552   cout << " Size of MRPC (cm) along Z                "<< AliTOFGeometry::ZPad()<<endl;
 553   cout << " Module Lengths (cm)" << endl;
 554   cout << " A Module: "<< AliTOFGeometry::ZlenA()<< "  B Modules: "<< AliTOFGeometry::ZlenB()<< "  C Modules: "<< AliTOFGeometry::ZlenC()<< endl;
 555   cout << " Inner radius of the TOF detector (cm): "<<AliTOFGeometry::Rmin() << endl;
 556   cout << " Outer radius of the TOF detector (cm): "<<AliTOFGeometry::Rmax() << endl;
 557   cout << " Max. half z-size of TOF (cm)         : "<<AliTOFGeometry::MaxhZtof() << endl;
 558   cout << " TOF Pad parameters   " << endl;
 559   cout << " Time Resolution (ns) "<< fTimeResolution <<" Pad Efficiency: "<< fpadefficiency << endl;
 560   cout << " Edge Effect option:  "<<  fEdgeEffect<< endl;
 561
 562   cout << " Boundary Effect Simulation Parameters " << endl;
 563   cout << " Hparameter: "<< fHparameter<<"  H2parameter:"<< fH2parameter <<"  Kparameter:"<< fKparameter<<"  K2parameter: "<< fK2parameter << endl;
 564   cout << " Efficiency in the central region of the pad: "<< fEffCenter << endl;
 565   cout << " Efficiency at the boundary region of the pad: "<< fEffBoundary << endl;
 566   cout << " Efficiency value at H2parameter "<< fEff2Boundary << endl;
 567   cout << " Efficiency value at K2parameter "<< fEff3Boundary << endl;
 568   cout << " Resolution (ps) in the central region of the pad: "<< fResCenter << endl;
 569   cout << " Resolution (ps) at the boundary of the pad      : "<< fResBoundary << endl;
 570   cout << " Slope (ps/K) for neighbouring pad               : "<< fResSlope <<endl;
 571   cout << " Time walk (ps) in the central region of the pad : "<< fTimeWalkCenter << endl;
 572   cout << " Time walk (ps) at the boundary of the pad       : "<< fTimeWalkBoundary<< endl;
 573   cout << " Slope (ps/K) for neighbouring pad               : "<< fTimeWalkSlope<<endl;
 574   cout << " Pulse Heigth Simulation Parameters " << endl;
 575   cout << " Flag for delay due to the PulseHeightEffect: "<< fTimeDelayFlag <<endl;
 576   cout << " Pulse Height Slope                           : "<< fPulseHeightSlope<<endl;
 577   cout << " Time Delay Slope                             : "<< fTimeDelaySlope<<endl;
 578   cout << " Minimum charge amount which could be induced : "<< fMinimumCharge<<endl;
 579   cout << " Smearing in charge in (q1/q2) vs x plot      : "<< fChargeSmearing<<endl;
 580   cout << " Smearing in log of charge ratio              : "<< fLogChargeSmearing<<endl;
 581   cout << " Smearing in time in time vs log(q1/q2) plot  : "<< fTimeSmearing<<endl;
 582   cout << " Flag for average time                        : "<< fAverageTimeFlag<<endl;
 583   cout << " Charge factor flag for matching              : "<< fChargeFactorForMatching<<endl;
 584   cout << " Edge tails option                            : "<< fEdgeTails << endl;
 585   cout << " TPC tracking  parameters " << endl;
 586   cout << " TPC tracking efficiency                      : "<< fTrackingEfficiency<< endl;
 587   cout << " Sigma vs momentum dependency flag            : "<< fSigmavsp << endl;
 588   cout << " Space uncertainties (cm). sigma(z) (cm): "<< fSigmaZ << " sigma(R(phi)) (cm): "<< fSigmarphi << endl;
 589   cout << " Momentum uncertainties.   sigma(delta(P)/P): "<< fSigmap <<" sigma(phi) (rad): "<< fSigmaPhi <<" sigma(theta) (rad): "<< fSigmaTheta << endl;
 590   cout << " Parameters for additional noise hits " << endl;
 591   cout << " Number of noise hits : " << fNoise <<" Slope parameter (ns) in the time distribution: " << fNoiseSlope << endl;
 592   cout << " Mean TOF for noise from outer regions (ns)" <<  fNoiseMeanTof << endl;
 593   cout << " Physical parameters " << endl;
 594   cout << " Magnetic Field (tesla)                   : "<< fField <<endl;
 595   cout << " Radiation length of the outer wall of TPC: "<< fRadLenTPC << endl;
 596   cout << " (TPC tracks)-(TOF pads) matching parameters " << endl;
 597   cout << " TRD Correction flag       : "<< fCorrectionTRD <<endl;
 598   cout << " Number of the last TPC row: "<< fLastTPCRow <<" Vertex radius (cm) for selected tracks: "<<fRadiusvtxBound<<endl;
 599   cout << " Max. number of test tracks: "<<fMaxTestTracks << endl;
 600   cout << " Space step (cm)           : "<< fStep <<endl;
 601   cout << " Matching style option     : "<< fMatchingStyle <<endl;
 602   cout << " Array parameters " << endl;
 603   cout << " Max.number of pads involved in the matching procedure: "<< fMaxPixels << endl;
 604   cout << " Max.number of TOF hits per event                     : "<< fMaxTOFHits<< endl;
 605   cout << " Max.number of tracks selected for matching           : "<< fMaxTracks << endl;
 606   cout << " Max.number of all tracks including the neutral ones  : "<< fMaxAllTracks<< endl;
 607   cout << " Debug Flag                                           : "<< fdbg << endl;
 608   cout << " Cut on momentum for selecting tracks                 : "<< fPBound << endl;
 609
 610 }
 611
 612 //__________________________________________________________________
 613 void AliTOFReconstructioner::IsInsideThePad(TVirtualMC *vmc, Float_t x, Float_t y, Float_t z, Int_t *nGeom, Float_t& zPad, Float_t& xPad)
 614 {
 615   //   input: x,y,z - coordinates of a hit
 616   //   output: array  nGeom[]
 617   //          nGeom[0] - the TOF sector number, 1,2,...,18 along azimuthal direction starting from -90 deg.!!!
 618   //          nGeom[1] - the TOF module number, 1,2,3,4,5=C,B,A,B,C along z-direction
 619   //          nGeom[2] - the TOF strip  number, 1,2,... along z-direction
 620   //          nGeom[3] - the TOF padz  number,  1,2=NPZ across a strip
 621   //          nGeom[4] - the TOF padx  number,  1,2,...,48=NPX along a strip
 622   //          zPad, xPad - coordinates of the hit in the pad frame
 623   //  numbering is adopted for the version 3.05 of AliRoot
 624   //  example:
 625   //   from Hits: sec,pla,str,padz,padx=4,2,14,2,35
 626   //  Vol. n.0: ALIC, copy number 1
 627   //  Vol. n.1: B077, copy number 1
 628   //  Vol. n.2: B074, copy number 5
 629   //  Vol. n.3: BTO2, copy number 1
 630   //  Vol. n.4: FTOB, copy number 2
 631   //  Vol. n.5: FLTB, copy number 0
 632   //  Vol. n.6: FSTR, copy number 14
 633   //  Vol. n.7: FSEN, copy number 0
 634   //  Vol. n.8: FSEZ, copy number 2
 635   //  Vol. n.9: FSEX, copy number 35
 636   //  Vol. n.10: FPAD, copy number 0
 637
 638
 639   Float_t xTOF[3];
 640   Int_t sector=0,module=0,strip=0,padz=0,padx=0;
 641   Int_t i,numed,nLevel,copyNumber;
 642   Gcvolu_t* gcvolu;
 643   char name[5];
 644   name[4]=0;
 645
 646   for (i=0; i<AliTOFGeometry::MaxTOFTree(); i++) nGeom[i]=0;
 647   zPad=100.;
 648   xPad=100.;
 649
 650   xTOF[0]=x;
 651   xTOF[1]=y;
 652   xTOF[2]=z;
 653
 654   TGeant3 * g3 = (TGeant3*) vmc;
 655
 656   g3->Gmedia(xTOF, numed);
 657   gcvolu=g3->Gcvolu();
 658   nLevel=gcvolu->nlevel;
 659   if(fdbg) {
 660     for (Int_t i=0; i<nLevel; i++) {
 661       strncpy(name,(char*) (&gcvolu->names[i]),4);
 662       cout<<"Vol. n."<<i<<": "<<name<<", copy number "<<gcvolu->number[i]<<endl;
 663     }
 664   }
 665   if(nLevel>=2) {
 666     // sector type name: B071(1,2,...,10),B074(1,2,3,4,5-PHOS),B075(1,2,3-RICH)
 667     strncpy(name,(char*) (&gcvolu->names[2]),4);
 668     // volume copy: 1,2,...,10 for B071, 1,2,3,4,5 for B074, 1,2,3 for B075
 669     copyNumber=gcvolu->number[2];
 670    if(!strcmp(name,"B071")) {
 671      if (copyNumber>=6 && copyNumber<=8) {
 672        sector=copyNumber+10;
 673      } else if (copyNumber>=1 && copyNumber<=5){
 674        sector=copyNumber+7;
 675      } else {
 676        sector=copyNumber-8;
 677      }
 678    } else if(!strcmp(name,"B075")) {
 679      sector=copyNumber+12;
 680    } else if(!strcmp(name,"B074")) {
 681      if (copyNumber>=1 && copyNumber<=3){
 682        sector=copyNumber+4;
 683      } else {
 684        sector=copyNumber-1;
 685      }
 686    }
 687   }
 688   if(sector) {
 689     nGeom[0]=sector;
 690     if(nLevel>=4) {
 691       // we'll use the module value in z-direction:
 692       //                                    1    2    3    4    5
 693       // the module order in z-direction: FTOC,FTOB,FTOA,FTOB,FTOC
 694       // the module copy:                   2    2    0    1    1
 695       // module type name: FTOA, FTOB, FTOC
 696       strncpy(name,(char*) (&gcvolu->names[4]),4);
 697       // module copy:
 698       copyNumber=gcvolu->number[4];
 699       if(!strcmp(name,"FTOC")) {
 700         if (copyNumber==2) {
 701           module=1;
 702         } else {
 703           module=5;
 704         }
 705       } else if(!strcmp(name,"FTOB")) {
 706         if (copyNumber==2) {
 707           module=2;
 708         } else {
 709           module=4;
 710         }
 711       } else if(!strcmp(name,"FTOA")) {
 712         module=3;
 713       }
 714     }
 715   }
 716
 717   if(module) {
 718     nGeom[1]=module;
 719     if(nLevel>=6) {
 720       // strip type name: FSTR
 721       strncpy(name,(char*) (&gcvolu->names[6]),4);
 722       // strip copy:
 723       copyNumber=gcvolu->number[6];
 724       if(!strcmp(name,"FSTR")) strip=copyNumber;
 725     }
 726   }
 727
 728   if(strip) {
 729     nGeom[2]=strip;
 730     if(nLevel>=8) {
 731       // padz type name: FSEZ
 732       strncpy(name,(char*) (&gcvolu->names[8]),4);
 733       // padz copy:
 734       copyNumber=gcvolu->number[8];
 735       if(!strcmp(name,"FSEZ")) padz=copyNumber;
 736     }
 737   }
 738   if(padz) {
 739     nGeom[3]=padz;
 740     if(nLevel>=9) {
 741       // padx type name: FSEX
 742       strncpy(name,(char*) (&gcvolu->names[9]),4);
 743       // padx copy:
 744       copyNumber=gcvolu->number[9];
 745       if(!strcmp(name,"FSEX")) padx=copyNumber;
 746     }
 747   }
 748
 749   if(padx) {
 750     nGeom[4]=padx;
 751     zPad=gcvolu->glx[2];  // check here
 752     xPad=gcvolu->glx[0];  // check here
 753   }
 754
 755   //   printf(" nGeom[0,1,2,3,4]=%i,%i,%i,%i,%i\n",nGeom[0],nGeom[1],nGeom[2],nGeom[3],nGeom[4]);
 756 }
 757
 758 //__________________________________________________________________
 759 void AliTOFReconstructioner::EpMulScatt(Float_t& px, Float_t& py, Float_t& pz, Float_t& p, Float_t& theta)
 760 {
 761   //   Momentum p  - before mult.scat.
 762   //   Momentum p2 - after mult.scat.
 763   //   THE0 - r.m.s. of deviation angle in plane
 764   //           (see RPP'96: Phys.Rev.D54 (1996) 134)
 765
 766   Float_t pt,thex,they,tantx,tanty,p2px,p2py,p2pz,costhe,sinthe,cospsi,sinpsi,p2x,p2y,p2z,p2,g;
 767
 768   pt=TMath::Sqrt(px*px+py*py);
 769   //   angles for p in the ' frame with Z'along p
 770   if(fMatchingStyle==1) {
 771     thex=theta*gRandom->Gaus();
 772     they=theta*gRandom->Gaus();
 773   } else {
 774     thex=3*(-theta+2*theta*gRandom->Rndm());
 775     they=3*(-theta+2*theta*gRandom->Rndm());
 776   }
 777   tantx=TMath::Tan(thex);
 778   tanty=TMath::Tan(they);
 779
 780   //   p2p - p2 in the ' frame
 781   p2pz=p/TMath::Sqrt(1.+tantx*tantx+tanty*tanty);
 782   p2py=p2pz*tanty;
 783   p2px=p2pz*tantx;
 784   //   choose X'so that PHI=0 (see Il'in, Pozdnyak Analiticheskaya geometriya, 1968, c.88
 785   //   for Euler angles PSI, THETA (PHI=0)
 786   costhe=pz/p;
 787   sinthe=pt/p;
 788   cospsi=-py/pt;
 789   sinpsi=px/pt;
 790   //
 791   g=p2py*costhe-p2pz*sinthe;
 792   p2x=p2px*cospsi-g*sinpsi;
 793   p2y=p2px*sinpsi+g*cospsi;
 794   p2z=p2py*sinthe+p2pz*costhe;
 795   p2=TMath::Sqrt(p2x*p2x+p2y*p2y+p2z*p2z);
 796
 797   //   Test angle
 798   g=(px*p2x+py*p2y+pz*p2z)/(p*p2);
 799   if(g>1) g=1;
 800   theta=TMath::ACos(g);
 801   px=p2x;
 802   py=p2y;
 803   pz=p2z;
 804   p=p2;
 805
 806 }
 807
 808 // std border effect algorithm
 809 //__________________________________________________________________
 810 void AliTOFReconstructioner::BorderEffect(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
 811 {
 812   // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
 813   //         geantTime - time generated by Geant, ns
 814   // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
 815   //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
 816   //         qInduced[iPad]- charge induced on pad, arb. units
 817   //                         this array is initialized at zero by the caller
 818   //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
 819   //                                   this array is initialized at zero by the caller
 820   //         averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
 821   //                       The weight is given by the qInduced[iPad]/qCenterPad
 822   //                                   this variable is initialized at zero by the caller
 823   //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
 824   //                                   this variable is initialized at zero by the caller
 825   //
 826   // Description of used variables:
 827   //         eff[iPad] - efficiency of the pad
 828   //         res[iPad] - resolution of the pad, ns
 829   //         timeWalk[iPad] - time walk of the pad, ns
 830   //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
 831   //         PadId[iPad] - Pad Identifier
 832   //                    E | F    -->   PadId[iPad] = 5 | 6
 833   //                    A | B    -->   PadId[iPad] = 1 | 2
 834   //                    C | D    -->   PadId[iPad] = 3 | 4
 835   //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
 836   //         qCenterPad - charge extimated for each pad, arb. units
 837   //         weightsSum - sum of weights extimated for each pad fired, arb. units
 838
 839   const Float_t kSigmaForTail[2] = {AliTOFGeometry::SigmaForTail1(),AliTOFGeometry::SigmaForTail2()}; //for tail
 840   Int_t iz = 0, ix = 0;
 841   Float_t dX = 0., dZ = 0., x = 0., z = 0.;
 842   Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
 843   Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
 844   Float_t logOfqInd = 0.;
 845   Float_t weightsSum = 0.;
 846   Int_t nTail[4]  = {0,0,0,0};
 847   Int_t padId[4]  = {0,0,0,0};
 848   Float_t eff[4]  = {0.,0.,0.,0.};
 849   Float_t res[4]  = {0.,0.,0.,0.};
 850   //  Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
 851   Float_t qCenterPad = 1.;
 852   Float_t timeWalk[4]  = {0.,0.,0.,0.};
 853   Float_t timeDelay[4] = {0.,0.,0.,0.};
 854
 855   nActivatedPads = 0;
 856   nFiredPads = 0;
 857
 858   (z0 <= 0) ? iz = 0 : iz = 1;
 859   dZ = z0 + (0.5 * AliTOFGeometry::NpadZ() - iz - 0.5) * AliTOFGeometry::ZPad(); // hit position in the pad frame, (0,0) - center of the pad
 860   z = 0.5 * AliTOFGeometry::ZPad() - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
 861   iz++;                                                                              // z row: 1, ..., AliTOFGeometry::NpadZ() = 2
 862   ix = (Int_t)((x0 + 0.5 * AliTOFGeometry::NpadX() * AliTOFGeometry::XPad()) / AliTOFGeometry::XPad());
 863   dX = x0 + (0.5 * AliTOFGeometry::NpadX() - ix - 0.5) * AliTOFGeometry::XPad(); // hit position in the pad frame, (0,0) - center of the pad
 864   x = 0.5 * AliTOFGeometry::XPad() - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
 865   ix++;                                                                              // x row: 1, ..., AliTOFGeometry::NpadX() = 48
 866
 867   ////// Pad A:
 868   nActivatedPads++;
 869   nPlace[nActivatedPads-1] = (iz - 1) * AliTOFGeometry::NpadX() + ix;
 870   qInduced[nActivatedPads-1] = qCenterPad;
 871   padId[nActivatedPads-1] = 1;
 872
 873   if (fEdgeEffect == 0) {
 874     eff[nActivatedPads-1] = fEffCenter;
 875     if (gRandom->Rndm() < eff[nActivatedPads-1]) {
 876       nFiredPads = 1;
 877       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns;
 878       isFired[nActivatedPads-1] = kTRUE;
 879       tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
 880       averageTime = tofTime[nActivatedPads-1];
 881     }
 882   } else {
 883
 884     if(z < h) {
 885       if(z < h2) {
 886         effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
 887       } else {
 888         effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
 889       }
 890       resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
 891       timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
 892       nTail[nActivatedPads-1] = 1;
 893     } else {
 894       effZ = fEffCenter;
 895       resZ = fResCenter;
 896       timeWalkZ = fTimeWalkCenter;
 897     }
 898
 899     if(x < h) {
 900       if(x < h2) {
 901         effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
 902       } else {
 903         effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
 904       }
 905       resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
 906       timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
 907       nTail[nActivatedPads-1] = 1;
 908     } else {
 909       effX = fEffCenter;
 910       resX = fResCenter;
 911       timeWalkX = fTimeWalkCenter;
 912     }
 913
 914     (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
 915     (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns
 916     (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 917
 918
 919     ////// Pad B:
 920     if(z < k2) {
 921       effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
 922     } else {
 923       effZ = fEff3Boundary * (k - z) / (k - k2);
 924     }
 925     resZ = fResBoundary + fResSlope * z / k;
 926     timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
 927
 928     if(z < k && z > 0) {
 929       if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
 930         nActivatedPads++;
 931         nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX();
 932         eff[nActivatedPads-1] = effZ;
 933         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
 934         timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
 935         nTail[nActivatedPads-1] = 2;
 936         if (fTimeDelayFlag) {
 937           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
 938           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
 939           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
 940           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
 941           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 942         } else {
 943           timeDelay[nActivatedPads-1] = 0.;
 944         }
 945         padId[nActivatedPads-1] = 2;
 946       }
 947     }
 948
 949
 950     ////// Pad C, D, E, F:
 951     if(x < k2) {
 952       effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
 953     } else {
 954       effX = fEff3Boundary * (k - x) / (k - k2);
 955     }
 956     resX = fResBoundary + fResSlope*x/k;
 957     timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
 958
 959     if(x < k && x > 0) {
 960       //   C:
 961       if(ix > 1 && dX < 0) {
 962         nActivatedPads++;
 963         nPlace[nActivatedPads-1] = nPlace[0] - 1;
 964         eff[nActivatedPads-1] = effX;
 965         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
 966         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 967         nTail[nActivatedPads-1] = 2;
 968         if (fTimeDelayFlag) {
 969           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
 970           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
 971           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
 972           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
 973           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 974         } else {
 975           timeDelay[nActivatedPads-1] = 0.;
 976         }
 977         padId[nActivatedPads-1] = 3;
 978
 979         //     D:
 980         if(z < k && z > 0) {
 981           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
 982             nActivatedPads++;
 983             nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() - 1;
 984             eff[nActivatedPads-1] = effX * effZ;
 985             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
 986             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 987
 988             nTail[nActivatedPads-1] = 2;
 989             if (fTimeDelayFlag) {
 990               if (TMath::Abs(x) < TMath::Abs(z)) {
 991                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
 992                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
 993                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
 994                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
 995               } else {
 996                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
 997                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
 998                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
 999                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1000               }
1001               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1002             } else {
1003               timeDelay[nActivatedPads-1] = 0.;
1004             }
1005             padId[nActivatedPads-1] = 4;
1006           }
1007         }  // end D
1008       }  // end C
1009
1010       //   E:
1011       if(ix < AliTOFGeometry::NpadX() && dX > 0) {
1012         nActivatedPads++;
1013         nPlace[nActivatedPads-1] = nPlace[0] + 1;
1014         eff[nActivatedPads-1] = effX;
1015         res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
1016         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1017         nTail[nActivatedPads-1] = 2;
1018         if (fTimeDelayFlag) {
1019           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1020           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1021           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
1022           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1023           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1024         } else {
1025           timeDelay[nActivatedPads-1] = 0.;
1026         }
1027         padId[nActivatedPads-1] = 5;
1028
1029
1030         //     F:
1031         if(z < k && z > 0) {
1032           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1033             nActivatedPads++;
1034             nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() + 1;
1035             eff[nActivatedPads - 1] = effX * effZ;
1036             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1037             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
1038             nTail[nActivatedPads-1] = 2;
1039             if (fTimeDelayFlag) {
1040               if (TMath::Abs(x) < TMath::Abs(z)) {
1041                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1042                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1043                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
1044                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1045               } else {
1046                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1047                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1048                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
1049                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1050               }
1051               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1052             } else {
1053               timeDelay[nActivatedPads-1] = 0.;
1054             }
1055             padId[nActivatedPads-1] = 6;
1056           }
1057         }  // end F
1058       }  // end E
1059     } // end if(x < k)
1060
1061
1062     for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
1063       if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
1064       if(gRandom->Rndm() < eff[iPad]) {
1065         isFired[iPad] = kTRUE;
1066         nFiredPads++;
1067         if(fEdgeTails) {
1068           if(nTail[iPad] == 0) {
1069             tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1070           } else {
1071             ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
1072             Double_t timeAB = ftail->GetRandom();
1073             tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
1074           }
1075         } else {
1076           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1077         }
1078         if (fAverageTimeFlag) {
1079           averageTime += tofTime[iPad] * qInduced[iPad];
1080           weightsSum += qInduced[iPad];
1081         } else {
1082           averageTime += tofTime[iPad];
1083           weightsSum += 1.;
1084         }
1085       }
1086     }
1087     if (weightsSum!=0) averageTime /= weightsSum;
1088   } // end else (fEdgeEffect != 0)
1089 }
1090
1091
1092 /* new algorithm (to be checked)
1093 //__________________________________________________________________
1094 void AliTOFReconstructioner::BorderEffect(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
1095 {
1096   // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
1097   //         geantTime - time generated by Geant, ns
1098   // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
1099   //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
1100   //         qInduced[iPad]- charge induced on pad, arb. units
1101   //                         this array is initialized at zero by the caller
1102   //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
1103   //                                   this array is initialized at zero by the caller
1104   //         averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
1105   //                       The weight is given by the qInduced[iPad]/qCenterPad
1106   //                                   this variable is initialized at zero by the caller
1107   //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
1108   //                                   this variable is initialized at zero by the caller
1109   //
1110   // Description of used variables:
1111   //         eff[iPad] - efficiency of the pad
1112   //         res[iPad] - resolution of the pad, ns
1113   //         timeWalk[iPad] - time walk of the pad, ns
1114   //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
1115   //         PadId[iPad] - Pad Identifier
1116   //                    E | F    -->   PadId[iPad] = 5 | 6
1117   //                    A | B    -->   PadId[iPad] = 1 | 2
1118   //                    C | D    -->   PadId[iPad] = 3 | 4
1119   //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
1120   //         qCenterPad - charge extimated for each pad, arb. units
1121   //         weightsSum - sum of weights extimated for each pad fired, arb. units
1122
1123   const Float_t kSigmaForTail[2] = {AliTOFGeometry::SigmaForTail1(),AliTOFGeometry::SigmaForTail2()}; //for tail
1124   Int_t iz = 0, ix = 0;
1125   Float_t dX = 0., dZ = 0., x = 0., z = 0.;
1126   Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
1127   Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
1128   Float_t logOfqInd = 0.;
1129   Float_t weightsSum = 0.;
1130   Int_t nTail[4]  = {0,0,0,0};
1131   Int_t padId[4]  = {0,0,0,0};
1132   Float_t eff[4]  = {0.,0.,0.,0.};
1133   Float_t res[4]  = {0.,0.,0.,0.};
1134   Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
1135   Float_t timeWalk[4]  = {0.,0.,0.,0.};
1136   Float_t timeDelay[4] = {0.,0.,0.,0.};
1137
1138   nActivatedPads = 0;
1139   nFiredPads = 0;
1140
1141   (z0 <= 0) ? iz = 0 : iz = 1;
1142   dZ = z0 + (0.5 * AliTOFGeometry::NpadZ() - iz - 0.5) * AliTOFGeometry::ZPad(); // hit position in the pad frame, (0,0) - center of the pad
1143   z = 0.5 * AliTOFGeometry::ZPad() - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
1144   iz++;                                                                              // z row: 1, ..., AliTOFGeometry::NpadZ() = 2
1145   ix = (Int_t)((x0 + 0.5 * AliTOFGeometry::NpadX() * AliTOFGeometry::XPad()) / AliTOFGeometry::XPad());
1146   dX = x0 + (0.5 * AliTOFGeometry::NpadX() - ix - 0.5) * AliTOFGeometry::XPad(); // hit position in the pad frame, (0,0) - center of the pad
1147   x = 0.5 * AliTOFGeometry::XPad() - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
1148   ix++;                                                                              // x row: 1, ..., AliTOFGeometry::NpadX() = 48
1149
1150   ////// Pad A:
1151   nActivatedPads++;
1152   nPlace[nActivatedPads-1] = (iz - 1) * AliTOFGeometry::NpadX() + ix;
1153   qInduced[nActivatedPads-1] = qCenterPad;
1154   padId[nActivatedPads-1] = 1;
1155
1156   if (fEdgeEffect == 0) {
1157     eff[nActivatedPads-1] = fEffCenter;
1158     if (gRandom->Rndm() < eff[nActivatedPads-1]) {
1159       nFiredPads = 1;
1160       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns;
1161       isFired[nActivatedPads-1] = kTRUE;
1162       tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
1163       averageTime = tofTime[nActivatedPads-1];
1164     }
1165   } else {
1166
1167     if(z < h) {
1168       if(z < h2) {
1169         effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
1170       } else {
1171         effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
1172       }
1173       resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
1174       timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
1175       nTail[nActivatedPads-1] = 1;
1176     } else {
1177       effZ = fEffCenter;
1178       resZ = fResCenter;
1179       timeWalkZ = fTimeWalkCenter;
1180     }
1181
1182     if(x < h) {
1183       if(x < h2) {
1184         effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
1185       } else {
1186         effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
1187       }
1188       resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
1189       timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
1190       nTail[nActivatedPads-1] = 1;
1191     } else {
1192       effX = fEffCenter;
1193       resX = fResCenter;
1194       timeWalkX = fTimeWalkCenter;
1195     }
1196
1197     (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
1198     (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2  ns
1199     (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1200
1201
1202     ////// Pad B:
1203     if(z < k2) {
1204       effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
1205     } else {
1206       effZ = fEff3Boundary * (k - z) / (k - k2);
1207     }
1208     resZ = fResBoundary + fResSlope * z / k;
1209     timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
1210
1211     if(z < k && z > 0) {
1212       if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1213         nActivatedPads++;
1214         nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX();
1215         eff[nActivatedPads-1] = effZ;
1216         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1217         timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
1218         nTail[nActivatedPads-1] = 2;
1219         if (fTimeDelayFlag) {
1220           qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1221           qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1222           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1223           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1224         } else {
1225           timeDelay[nActivatedPads-1] = 0.;
1226         }
1227         padId[nActivatedPads-1] = 2;
1228       }
1229     }
1230
1231
1232     ////// Pad C, D, E, F:
1233     if(x < k2) {
1234       effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
1235     } else {
1236       effX = fEff3Boundary * (k - x) / (k - k2);
1237     }
1238     resX = fResBoundary + fResSlope*x/k;
1239     timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
1240
1241     if(x < k && x > 0) {
1242       //   C:
1243       if(ix > 1 && dX < 0) {
1244         nActivatedPads++;
1245         nPlace[nActivatedPads-1] = nPlace[0] - 1;
1246         eff[nActivatedPads-1] = effX;
1247         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1248         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1249         nTail[nActivatedPads-1] = 2;
1250         if (fTimeDelayFlag) {
1251           qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1252           qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1253           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1254           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1255         } else {
1256           timeDelay[nActivatedPads-1] = 0.;
1257         }
1258         padId[nActivatedPads-1] = 3;
1259
1260         //     D:
1261         if(z < k && z > 0) {
1262           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1263             nActivatedPads++;
1264             nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() - 1;
1265             eff[nActivatedPads-1] = effX * effZ;
1266             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1267             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1268
1269             nTail[nActivatedPads-1] = 2;
1270             if (fTimeDelayFlag) {
1271               if (TMath::Abs(x) < TMath::Abs(z)) {
1272                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1273                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1274                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1275               } else {
1276                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1277                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1278                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1279               }
1280               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1281             } else {
1282               timeDelay[nActivatedPads-1] = 0.;
1283             }
1284             padId[nActivatedPads-1] = 4;
1285           }
1286         }  // end D
1287       }  // end C
1288
1289       //   E:
1290       if(ix < AliTOFGeometry::NpadX() && dX > 0) {
1291         nActivatedPads++;
1292         nPlace[nActivatedPads-1] = nPlace[0] + 1;
1293         eff[nActivatedPads-1] = effX;
1294         res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
1295         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1296         nTail[nActivatedPads-1] = 2;
1297         if (fTimeDelayFlag) {
1298           qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1299           qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1300           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1301           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1302         } else {
1303           timeDelay[nActivatedPads-1] = 0.;
1304         }
1305         padId[nActivatedPads-1] = 5;
1306
1307
1308         //     F:
1309         if(z < k && z > 0) {
1310           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1311             nActivatedPads++;
1312             nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() + 1;
1313             eff[nActivatedPads - 1] = effX * effZ;
1314             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1315             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
1316             nTail[nActivatedPads-1] = 2;
1317             if (fTimeDelayFlag) {
1318               if (TMath::Abs(x) < TMath::Abs(z)) {
1319                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1320                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1321                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1322               } else {
1323                 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1324                 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1325                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1326               }
1327               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1328             } else {
1329               timeDelay[nActivatedPads-1] = 0.;
1330             }
1331             padId[nActivatedPads-1] = 6;
1332           }
1333         }  // end F
1334       }  // end E
1335     } // end if(x < k)
1336
1337
1338     for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
1339       if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
1340       if(gRandom->Rndm() < eff[iPad]) {
1341         isFired[iPad] = kTRUE;
1342         nFiredPads++;
1343         if(fEdgeTails) {
1344           if(nTail[iPad] == 0) {
1345             tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1346           } else {
1347             ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
1348             Double_t timeAB = ftail->GetRandom();
1349             tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
1350           }
1351         } else {
1352           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1353         }
1354         if (fAverageTimeFlag) {
1355           averageTime += tofTime[iPad] * qInduced[iPad];
1356           weightsSum += qInduced[iPad];
1357         } else {
1358           averageTime += tofTime[iPad];
1359           weightsSum += 1.;
1360         }
1361       }
1362     }
1363     if (weightsSum!=0) averageTime /= weightsSum;
1364
1365   } // end else (fEdgeEffect != 0)
1366
1367   //cout << "timedelay " << timeDelay[0] << endl;
1368   //cout << "timedelay " << timeDelay[1] << endl;
1369   //cout << "timedelay " << timeDelay[2] << endl;
1370   //cout << "timedelay " << timeDelay[3] << endl;
1371
1372 }
1373 */
1374
1375
1376 //__________________________________________________________________
1377 Int_t AliTOFReconstructioner::PDGtoGeantCode(Int_t pdgcode)
1378 {
1379   //
1380   // Gives the GEANT code from KF code of LUND JETSET
1381   //
1382   Int_t geantCode=0; // default value
1383   switch (pdgcode) {
1384   case 22:
1385     geantCode=1;            // GAMMA
1386     break ;
1387   case -11:
1388     geantCode=2;            // E+
1389     break ;
1390   case 11:
1391     geantCode=3;            // E-
1392     break ;
1393   case 12:
1394     geantCode=4;            // NUE
1395     break ;
1396   case 14:
1397     geantCode=4;            // NUMU
1398     break ;
1399   case -13:
1400     geantCode=5;            // MU+
1401     break ;
1402   case 13:
1403     geantCode=6;            // MU-
1404     break ;
1405   case 111:
1406     geantCode=7;            // PI0
1407     break ;
1408   case 211:
1409     geantCode=8;            // PI+
1410     break ;
1411   case -211:
1412     geantCode=9;            // PI-
1413     break ;
1414   case 130:
1415     geantCode=10;           // K_L0
1416     break ;
1417   case 321:
1418     geantCode=11;           // K+
1419     break ;
1420   case -321:
1421     geantCode=12;           // K-
1422     break ;
1423   case 2112:
1424     geantCode=13;           // N0
1425     break ;
1426   case 2212:
1427     geantCode=14;           // P+
1428     break ;
1429   case -2212:
1430     geantCode=15;           // P~-
1431     break ;
1432   case 310:
1433     geantCode=16;           // K_S0
1434     break ;
1435   case 221:
1436     geantCode=17;           // ETA
1437     break ;
1438   case 3122:
1439     geantCode=18;           // LAMBDA0
1440     break ;
1441   case 3222:
1442     geantCode=19;           // SIGMA+
1443     break ;
1444   case 3212:
1445     geantCode=20;           // SIGMA0
1446     break ;
1447   case 3112:
1448     geantCode=21;           // SIGMA-
1449     break ;
1450   case 3322:
1451     geantCode=22;           // XI0
1452     break ;
1453   case 3312:
1454     geantCode=23;           // XI-
1455     break ;
1456   case 3334:
1457     geantCode=24;           // OMEGA-
1458     break ;
1459   case -2112:
1460     geantCode=25;           // N~0
1461     break ;
1462   case -3122:
1463     geantCode=26;           // LAMBDA~0
1464     break ;
1465   case -3112:
1466     geantCode=27;           // SIGMA~+
1467     break ;
1468   case -3212:
1469     geantCode=28;           // SIGMA~0
1470     break ;
1471   case -3222:
1472     geantCode=29;           // SIGMA~-
1473     break ;
1474   case -3322:
1475     geantCode=30;           // XI~0
1476     break ;
1477   case -3312:
1478     geantCode=31;           // XI~+
1479     break ;
1480   case -3334:
1481     geantCode=32;           // OMEGA~+
1482     break ;
1483   case 223:
1484     geantCode=33;           // OMEGA(782)
1485     break ;
1486   case 333:
1487     geantCode=34;           // PHI(1020)
1488     break ;
1489   case 411:
1490     geantCode=35;           // D+
1491     break ;
1492   case -411:
1493     geantCode=36;           // D-
1494     break ;
1495   case 421:
1496     geantCode=37;           // D0
1497     break ;
1498   case -421:
1499     geantCode=38;           // D~0
1500     break ;
1501   case 431:
1502     geantCode=39;           // D_S+
1503     break ;
1504   case -431:
1505     geantCode=40;           // D_S~-
1506     break ;
1507   case 4122:
1508     geantCode=41;           // LAMBDA_C+
1509     break ;
1510   case 213:
1511     geantCode=42;           // RHP(770)+
1512     break ;
1513   case -213:
1514     geantCode=43;           // RHO(770)-
1515     break ;
1516   case 113:
1517     geantCode=44;           // RHO(770)0
1518     break ;
1519   default:
1520     geantCode=45;
1521     break;
1522   }
1523
1524   return geantCode;
1525 }
1526
1527 //__________________________________________________________________
1528 Bool_t AliTOFReconstructioner::operator==( AliTOFReconstructioner const & tofrec)const
1529 {
1530   // Equal operator.
1531   // Reconstructioners are equal if their parameters are equal
1532
1533   // split the member variables in analogous categories
1534
1535   // time resolution and edge effect parameters
1536   Bool_t dummy0=(fTimeResolution==tofrec.fTimeResolution)&&
1537           (fpadefficiency==tofrec.fpadefficiency)&&
1538           (fEdgeEffect==tofrec.fEdgeEffect)&&
1539           (fEdgeTails==tofrec.fEdgeTails)&&
1540           (fHparameter==tofrec.fHparameter)&&
1541           (fH2parameter==tofrec.fH2parameter)&&
1542           (fKparameter==tofrec.fKparameter)&&
1543           (fK2parameter==tofrec.fK2parameter);
1544
1545   // pad efficiency parameters
1546   Bool_t dummy1=(fEffCenter==tofrec.fEffCenter)&&
1547          (fEffBoundary==tofrec.fEffBoundary)&&
1548          (fEff2Boundary==tofrec.fEff2Boundary)&&
1549          (fEff3Boundary==tofrec.fEff3Boundary)&&
1550          (fResCenter==tofrec.fResCenter)&&
1551          (fResBoundary==tofrec.fResBoundary)&&
1552          (fResSlope==tofrec.fResSlope);
1553
1554   // time walk parameters
1555   Bool_t dummy2=(fTimeWalkCenter==tofrec.fTimeWalkCenter)&&
1556           (fTimeWalkBoundary==tofrec.fTimeWalkBoundary)&&
1557           (fTimeWalkSlope==tofrec.fTimeWalkSlope)&&
1558           (fTimeDelayFlag==tofrec.fTimeDelayFlag)&&
1559           (fPulseHeightSlope==tofrec.fPulseHeightSlope)&&
1560           (fTimeDelaySlope==tofrec.fTimeDelaySlope);
1561
1562   // ADC-TDC correlation parameters
1563   Bool_t dummy3=(fMinimumCharge==tofrec.fMinimumCharge)&&
1564           (fChargeSmearing==tofrec.fChargeSmearing )&&
1565           (fLogChargeSmearing==tofrec.fLogChargeSmearing )&&
1566           (fTimeSmearing==tofrec.fTimeSmearing )&&
1567           (fAverageTimeFlag==tofrec.fAverageTimeFlag)&&
1568           (fChargeFactorForMatching==tofrec.fChargeFactorForMatching)&&
1569           (fMatchingStyle==tofrec.fMatchingStyle);
1570
1571   Bool_t dummy4=(fTrackingEfficiency==tofrec.fTrackingEfficiency)&&
1572           (fSigmavsp==tofrec.fSigmavsp)&&
1573           (fSigmaZ==tofrec.fSigmaZ)&&
1574           (fSigmarphi==tofrec.fSigmarphi)&&
1575           (fSigmap==tofrec.fSigmap)&&
1576           (fSigmaPhi==tofrec.fSigmaPhi)&&
1577           (fSigmaTheta==tofrec.fSigmaTheta)&&
1578           (fNoise==tofrec.fNoise)&&
1579           (fNoiseSlope==tofrec.fNoiseSlope)&&
1580           (fField==tofrec.fField)&&
1581           (fRadLenTPC==tofrec.fRadLenTPC)&&
1582           (fCorrectionTRD==tofrec.fCorrectionTRD)&&
1583           (fLastTPCRow==tofrec.fLastTPCRow)&&
1584           (fRadiusvtxBound==tofrec.fRadiusvtxBound)&&
1585           (fMaxTestTracks==tofrec.fMaxTestTracks)&&
1586           (fStep==tofrec.fStep)&&
1587           (fMaxPixels==tofrec.fMaxPixels)&&
1588           (fMaxAllTracks==tofrec.fMaxAllTracks)&&
1589           (fMaxTracks==tofrec.fMaxTracks)&&
1590           (fMaxTOFHits==tofrec.fMaxTOFHits)&&
1591           (fPBound==tofrec.fPBound);
1592
1593   if( dummy0 && dummy1 && dummy2 && dummy3 && dummy4)
1594     return kTRUE ;
1595   else
1596     return kFALSE ;
1597
1598 }
1599 //____________________________________________________________________________
1600 void AliTOFReconstructioner::UseHitsFrom(const char * filename)
1601 {
1602   SetTitle(filename) ;
1603 }
1604
1605 //____________________________________________________________________________
1606 void AliTOFReconstructioner::InitArray(Float_t array[], Int_t nlocations)
1607 {
1608   //
1609   // Initialize the array of Float_t
1610   //
1611   for (Int_t i = 0; i < nlocations; i++) {
1612     array[i]=0.;
1613   }                             // end loop
1614
1615 }
1616
1617 //____________________________________________________________________________
1618 void AliTOFReconstructioner::InitArray(Int_t array[], Int_t nlocations)
1619 {
1620   //
1621   // Initialize the array of Int_t
1622   //
1623   for (Int_t i = 0; i < nlocations; i++) {
1624     array[i]=0;
1625   }                             // end loop
1626
1627 }
1628
1629
1630 //____________________________________________________________________________
1631 void AliTOFReconstructioner::ReadTOFHits(Int_t ntracks, TTree* treehits,
1632                 TClonesArray* tofhits, Int_t ***MapPixels, Int_t* kTOFhitFirst,
1633                 AliTOFPad* pixelArray , Int_t* iTOFpixel, Float_t* toftime,
1634                 AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel)
1635 {
1636   //
1637   // Read TOF hits for the current event and fill arrays
1638   //
1639   // Start loop on primary tracks in the hits containers
1640   //
1641   // Noise meaning in ReadTOFHits: we use the word 'noise' in the
1642   // following cases
1643   // - signals produced by secondary particles
1644   // - signals produced by the next hits (out of the first) of a given track
1645   //   (both primary and secondary)
1646   // - signals produced by edge effect
1647
1648
1649   TParticle *particle;
1650   Int_t nHitOutofTofVolumes; // number of hits out of TOF GEANT volumes (it happens in very
1651                              // few cases)
1652   Int_t * npixel = new Int_t[AliTOFGeometry::MaxTOFTree()]; // array used by TOFRecon for check on TOF geometry
1653   Int_t npions=0;    // number of pions for the current event
1654   Int_t nkaons=0;    // number of kaons for the current event
1655   Int_t nprotons=0;  // number of protons for the current event
1656   Int_t nelectrons=0;// number of electrons for the current event
1657   Int_t nmuons=0;    // number of muons for the current event
1658   Float_t tofpos[3];     // TOF hit position and GEANT time
1659   Float_t zPad,xPad;
1660   Int_t nbytes = 0;
1661   Int_t ipart, nhits=0, nHitsFromPrimaries=0;
1662   Int_t ntotalTOFhits=0; // total number of TOF hits for the current event
1663   Int_t ipartLast=-1;    // last track identifier
1664   Int_t iFirstHit;       // flag to check if the current hit is the first hit on TOF for the
1665                          // current track
1666   Int_t iNoiseHit=0;     // flag used to tag noise hits (the noise meaning is reported in the
1667                          // header of the ReadTOFHits method)
1668   Int_t nhitWithoutNoise;// number of hits not due to noise
1669   Int_t inoise=0,inoise2=0;
1670   Int_t nMultipleSignOnSamePad=0; // number of cases where a pad is fired more than one time
1671   Int_t nPixEdge=0;      // additional pads fired due to edge effect in ReadTOFHits (local var)
1672   // array used for counting different types of primary particles
1673   Int_t particleTypeGEANT[50]={0,4,4,0,5,5,0,3,3,0,
1674                    2,2,0,1,1,0,0,0,0,0,
1675                    0,0,0,0,0,0,0,0,0,0,
1676                    0,0,0,0,0,0,0,0,0,0,
1677                    0,0,0,0,0,0,0,0,0,0};
1678   Int_t particleType,particleInTOFtype[6][3];
1679   for (Int_t i=0;i<6;i++) {
1680     for (Int_t j=0;j<3;j++) {
1681       particleInTOFtype[i][j]=0;
1682     }
1683   }
1684
1685   // speed-up the code
1686   treehits->SetBranchStatus("*",0); // switch off all branches
1687   treehits->SetBranchStatus("TOF*",1); // switch on only TOF
1688
1689   for (Int_t track=0; track<ntracks;track++) { // starting loop on primary tracks for the current event
1690
1691     gAlice->ResetHits();
1692     nbytes += treehits->GetEvent(track);
1693     nhits = tofhits->GetEntriesFast();
1694
1695     ntotalTOFhits+=nhits;
1696
1697     // Start loop on hits connected to the current primary tracked particle
1698     // (including hits produced by secondary particles generaterd by the
1699     // current ptimary tracked particle)
1700     for (Int_t hit=0;hit<nhits;hit++) {
1701       AliTOFhit* tofHit   = (AliTOFhit*)tofhits->UncheckedAt(hit);
1702       ipart    = tofHit->GetTrack();
1703       if(ipart>=fMaxAllTracks) break;
1704       Float_t geantTime= tofHit->GetTof(); // it is given in [s]
1705       particle = (TParticle*)gAlice->GetMCApp()->Particle(ipart);
1706
1707       Int_t pdgCode=particle->GetPdgCode();
1708       // Only high momentum tracks (see fPBound value)
1709       // momentum components at vertex
1710       Float_t pxvtx = particle->Px();
1711       Float_t pyvtx = particle->Py();
1712       Float_t pzvtx = particle->Pz();
1713       Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx);
1714       if(pvtx>fPBound) {
1715
1716         if(particle->GetFirstMother() < 0) nHitsFromPrimaries++; // count primaries
1717
1718         // x and y coordinates of the particle production vertex
1719         Float_t vx = particle->Vx();
1720         Float_t vy = particle->Vy();
1721         Float_t vr = TMath::Sqrt(vx*vx+vy*vy); // cylindrical radius of the particle production vertex
1722
1723         Float_t x = tofHit->X(); tofpos[0]=x;
1724         Float_t y = tofHit->Y(); tofpos[1]=y;
1725         Float_t z = tofHit->Z(); tofpos[2]=z;
1726         /* var used for QA
1727         Float_t tofradius = TMath::Sqrt(x*x+y*y); // radius cilindrical coordinate of the TOF hit
1728         */
1729         // momentum components (cosine) when striking the TOF
1730         Float_t pxtof = tofHit->GetPx();
1731         Float_t pytof = tofHit->GetPy();
1732         Float_t pztof = tofHit->GetPz();
1733         // scalar product indicating the direction of the particle when striking the TOF
1734         /* var used for QA
1735         // (>0 for outgoing particles)
1736         Float_t isGoingOut = (x*pxtof+y*pytof+z*pztof)/TMath::Sqrt(x*x+y*y+z*z);
1737         */
1738         Float_t momtof = tofHit->GetMom();
1739         // now momentum components when striking the TOF
1740         pxtof *= momtof;
1741         pytof *= momtof;
1742         pztof *= momtof;
1743         particleType=particleTypeGEANT[PDGtoGeantCode(pdgCode)-1];
1744         if(particleType) {
1745           particleInTOFtype[5][2]++;
1746           particleInTOFtype[particleType-1][2]++;
1747         }
1748         iFirstHit=0;
1749         //  without noise hits
1750
1751         if(ipart!=ipartLast) {
1752           iFirstHit=1;
1753           toftime[ipart]=geantTime;    //time [s]
1754           // tofMom[ipart]=momtof;
1755           ipartLast=ipart;
1756           if(particle->GetFirstMother() < 0) {
1757             Int_t abspdgCode=TMath::Abs(pdgCode);
1758             switch (abspdgCode) {
1759             case 211:
1760               npions++;
1761               break ;
1762             case 321:
1763               nkaons++;
1764               break ;
1765             case 2212:
1766               nprotons++;
1767               break ;
1768             case 11:
1769               nelectrons++;
1770               break ;
1771             case 13:
1772               nmuons++;
1773               break ;
1774             }
1775           }
1776           if(vr>fRadiusvtxBound) {
1777             if(particleType) {
1778               particleInTOFtype[5][1]++;
1779               particleInTOFtype[particleType-1][1]++;
1780             }
1781             inoise++;
1782             inoise2++;
1783           } else {
1784             if(particleType) {
1785               particleInTOFtype[5][0]++;
1786               particleInTOFtype[particleType-1][0]++;
1787             }
1788           }
1789         } else {
1790           inoise++;
1791           if(particleType) {
1792             particleInTOFtype[5][1]++;
1793             particleInTOFtype[particleType-1][1]++;
1794           }
1795         } //end if(ipart!=ipartLast)
1796
1797         IsInsideThePad(gMC,x,y,z,npixel,zPad,xPad);
1798
1799         Int_t sec  = tofHit->GetSector();
1800         Int_t pla  = tofHit->GetPlate();
1801         Int_t str  = tofHit->GetStrip();
1802         if(sec!=npixel[0] || pla!=npixel[1] || str!=npixel[2]){// check on volume
1803           cout << "sector" << sec << " npixel[0] " << npixel[0] << endl;
1804           cout << "plate " << pla << " npixel[1] " << npixel[1] << endl;
1805           cout << "strip " << str << " npixel[2] " << npixel[2] << endl;
1806         } // close check on volume
1807
1808         Int_t padz = tofHit->GetPadz();
1809         Int_t padx = tofHit->GetPadx();
1810         Float_t Zpad = tofHit->GetDz();
1811         Float_t Xpad = tofHit->GetDx();
1812
1813
1814         if (npixel[4]==0){
1815           IsInsideThePad(gMC,x,y,z,npixel,zPad,xPad);
1816           if (npixel[4]==0){
1817             nHitOutofTofVolumes++;
1818           }
1819         } else {
1820           Float_t zStrip=AliTOFGeometry::ZPad()*(padz-0.5-0.5*AliTOFGeometry::NpadZ())+Zpad;
1821           if(padz!=npixel[3]) printf("            : Zpad=%f, padz=%i, npixel[3]=%i, zStrip=%f\n",Zpad,padz,npixel[3],zStrip);
1822           Float_t xStrip=AliTOFGeometry::XPad()*(padx-0.5-0.5*AliTOFGeometry::NpadX())+Xpad;
1823
1824           Int_t nPlace[4]={0,0,0,0};
1825           nPlace[0]=(padz-1)*AliTOFGeometry::NpadX()+padx;
1826
1827           Int_t   nActivatedPads=0;
1828           Int_t   nFiredPads=0;
1829           Bool_t  isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE};
1830           Float_t tofAfterSimul[4]={0.,0.,0.,0.};
1831           Float_t qInduced[4]={0.,0.,0.,0.};
1832           Float_t averageTime=0.;
1833
1834
1835           BorderEffect(zStrip,xStrip,geantTime*1.0e+09,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect
1836
1837
1838           if(nFiredPads) {
1839             for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
1840               if(isFired[indexOfPad]){// the pad has fired
1841                 if(indexOfPad==0) {// the hit belongs to a fired pad
1842                   isHitOnFiredPad++;
1843                   hitArray[isHitOnFiredPad-1].SetHit(ipart,pdgCode,tofpos,momtof,vr,iFirstHit);
1844                   iNoiseHit=0;
1845
1846                   if(vr>fRadiusvtxBound || iFirstHit==0) iNoiseHit=1;
1847
1848                   hitArray[isHitOnFiredPad-1].SetNoise(iNoiseHit);
1849                   if(iFirstHit) kTOFhitFirst[ipart]=isHitOnFiredPad;
1850
1851                 }// close - the hit belongs to a fired pad
1852
1853                 Int_t iMapFirstIndex=AliTOFGeometry::NSectors()*(npixel[1]-1)+npixel[0]-1;
1854                 Int_t iMapValue=MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1];
1855
1856                 if(iMapValue==0) {
1857                   ipixel++;
1858                   if(indexOfPad) {
1859                     iNoiseHit=1;
1860                     nPixEdge++;
1861                   } else {
1862                     iTOFpixel[ipart]=ipixel;
1863                   }
1864
1865                   if(ipixel>fMaxPixels){ // check on the total number of activated pads
1866                     cout << "ipixel=" << ipixel << " > fMaxPixels=" << fMaxPixels << endl;
1867                     return;
1868                   } // close check on the number of activated pads
1869
1870                   MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1]=ipixel;
1871                   pixelArray[ipixel-1].SetGeom(npixel[0],npixel[1],npixel[2],nPlace[indexOfPad]);
1872                   pixelArray[ipixel-1].SetTrack(ipart);
1873                   if(iNoiseHit) {
1874                     pixelArray[ipixel-1].AddState(1);
1875                   } else {
1876                     if(tofAfterSimul[indexOfPad]<0) cout << "Time of Flight after detector simulation is negative" << endl;
1877                     pixelArray[ipixel-1].AddState(10);
1878                   }
1879
1880                   pixelArray[ipixel-1].SetTofChargeHit(tofAfterSimul[indexOfPad],qInduced[indexOfPad],geantTime*1.0e+09,isHitOnFiredPad);
1881                 } else { //else if(iMapValue==0)
1882                   if(indexOfPad==0) iTOFpixel[ipart]=iMapValue;
1883                   nMultipleSignOnSamePad++;
1884
1885                   if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) {
1886                     pixelArray[iMapValue-1].SetTrack(ipart);
1887                     //                   if(indexOfPad==0) pixelArray[iMapValue-1].SetTrack(ipart);
1888                     if(indexOfPad) iNoiseHit=1;
1889                     if(iNoiseHit) {
1890                       pixelArray[iMapValue-1].AddState(1);
1891                     } else {
1892                       pixelArray[iMapValue-1].AddState(10);
1893                     }
1894                     pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]);
1895                     pixelArray[iMapValue-1].SetGeantTime(geantTime*1.0e+09);
1896                     pixelArray[iMapValue-1].SetHit(isHitOnFiredPad);
1897                   } // close if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetTime() )
1898                 }  //end of Pixel filling
1899               }  // close if(isFired[indexOfPad])
1900             }  //end loop on activated pads indexOfPad
1901           } // close if(nFiredPads)
1902         }  //end of hit with npixel[3]!=0
1903       }  //high momentum tracks
1904     }  //end on TOF hits
1905   }  //end on primary tracks
1906
1907
1908   if(fdbg) {
1909     cout << ntotalTOFhits << " - total number of TOF hits   " << nHitsFromPrimaries << " -  primary     " <<  endl;
1910     cout << inoise << " - noise hits, " << inoise2<< " - first crossing of a track with Rvtx>" << fRadiusvtxBound << endl;
1911     //   cout << inoise << " - noise hits (" << 100.*inoise/ihit << " %), " << inoise2
1912     //<< " - first crossing of a track with Rvtx>" << RVTXBOUND << endl;
1913     nhitWithoutNoise=isHitOnFiredPad;
1914
1915     cout << ipixel << " fired pixels (" << nMultipleSignOnSamePad << " multiple fired pads, " << endl;
1916     //j << " fired by noise, " << j1 << " noise+track)" <<  endl;
1917     printf(" %i additional pads are fired due to edge effect\n",nPixEdge);
1918     cout << npions <<   "   primary pions     reached TOF" << endl;
1919     cout << nkaons <<   "   primary kaons     reached TOF" << endl;
1920     cout << nprotons << "   primary protons   reached TOF" << endl;
1921     cout << nelectrons<<"   primary electrons reached TOF" << endl;
1922     cout << nmuons <<   "   primary muons     reached TOF" << endl;
1923     cout << "number of TOF hits for different species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl;
1924     cout << "   first number - track hits, second - noise ones, third - all" << endl;
1925     for (Int_t i=0;i<6;i++) cout << i+1 << "  " << particleInTOFtype[i][0] << "  " << particleInTOFtype[i][1] << "  " << particleInTOFtype[i][2] << endl;
1926
1927     Int_t primaryReachedTOF[6];
1928     primaryReachedTOF[0]=npions;
1929     primaryReachedTOF[1]=nkaons;
1930     primaryReachedTOF[2]=nprotons;
1931     primaryReachedTOF[3]=nelectrons;
1932     primaryReachedTOF[4]=nmuons;
1933     primaryReachedTOF[5]=npions+nkaons+nprotons+nelectrons+nmuons;
1934
1935     cout << " Reading TOF hits done" << endl;
1936   }
1937
1938   delete [] npixel;
1939 }
1940
1941 //____________________________________________________________________________
1942 void AliTOFReconstructioner::AddNoiseFromOuter(Option_t *option, Int_t ***MapPixels, AliTOFPad* pixelArray , AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel)
1943 {
1944   //
1945   // Add noise hits from outer regions (forward and backward) according
1946   // to parameterized fZNoise distribution (to be used with events
1947   // generated in the barrel region)
1948
1949   Float_t * zLen = new Float_t[AliTOFGeometry::NPlates()+1];
1950   Float_t * zStrips = new Float_t[AliTOFGeometry::NPlates()];
1951   zStrips[0]=(Float_t) (AliTOFGeometry::NStripC());
1952   zStrips[1]=(Float_t) (AliTOFGeometry::NStripB());
1953   zStrips[2]=(Float_t) (AliTOFGeometry::NStripA());
1954   zStrips[3]=(Float_t) (AliTOFGeometry::NStripB());
1955   zStrips[4]=(Float_t) (AliTOFGeometry::NStripC());
1956
1957   zLen[5]=AliTOFGeometry::ZlenA()*0.5+AliTOFGeometry::ZlenB()+AliTOFGeometry::ZlenC();
1958   zLen[4]=zLen[5]-AliTOFGeometry::ZlenC();
1959   zLen[3]=zLen[4]-AliTOFGeometry::ZlenB();
1960   zLen[2]=zLen[3]-AliTOFGeometry::ZlenA();
1961   zLen[1]=zLen[2]-AliTOFGeometry::ZlenB();
1962   zLen[0]=zLen[1]-AliTOFGeometry::ZlenC();
1963
1964
1965   Int_t isector;    // random sector number
1966   Int_t iplate;     // random plate number
1967   Int_t istrip;     // random strip number in the plate
1968   Int_t ipadAlongX; // random pad number along x direction
1969   Int_t ipadAlongZ; // random pad number along z direction
1970   Int_t ipad;
1971   Int_t nPixEdge=0; // additional pads fired due to edge effect when adding noise from outer
1972                     // regions
1973
1974   // x -> time of flight given in ns
1975   TF1 *noiseTof = new TF1("noiseTof","exp(-x/20)",0,100);
1976
1977   if(strstr(option,"pp")){
1978     fZnoise = new TF1("fZnoise","257.8-0.178*x-0.000457*x*x",-AliTOFGeometry::MaxhZtof(),AliTOFGeometry::MaxhZtof());
1979   }
1980   if(strstr(option,"Pb-Pb")){
1981     fZnoise = new TF1("fZnoise","182.2-0.09179*x-0.0001931*x*x",-AliTOFGeometry::MaxhZtof(),AliTOFGeometry::MaxhZtof());
1982   }
1983
1984   if(fNoise) {
1985     if(fdbg) cout << " Start adding additional noise  hits from outer regions" << endl;
1986
1987     for(Int_t i=0;i<fNoise;i++) {
1988
1989       isector=(Int_t) (AliTOFGeometry::NSectors()*gRandom->Rndm())+1; //the sector number
1990       //  non-flat z-distribution of additional hits
1991       Float_t zNoise=fZnoise->GetRandom();
1992
1993       // holes for PHOS and HMPID
1994       if(((AliTOF *) gAlice->GetDetector("TOF"))->IsVersion()==2) {
1995         // to be checked the holes case
1996         if(isector>12 && isector<16) { // sectors 13,14,15 - RICH
1997           do {
1998             iplate=(Int_t) (AliTOFGeometry::NPlates()*gRandom->Rndm())+1;
1999           } while (iplate==2 || iplate==3 || iplate==4);
2000           //         } else if(isector>11 && isector<17) { // sectors 12,13,14,15,16 - PHOS
2001         } else if(isector>2 && isector<8) { // sectors 3,4,5,6,7 - PHOS
2002           do {
2003             iplate=(Int_t) (AliTOFGeometry::NPlates()*gRandom->Rndm())+1;
2004           } while (iplate==3);
2005         } else {
2006           iplate=(Int_t) (AliTOFGeometry::NPlates()*gRandom->Rndm())+1;
2007         }
2008       } else {
2009         iplate=0;
2010         do {
2011           iplate++;
2012         } while(zNoise>zLen[iplate]);
2013       }
2014       // end of holes
2015
2016       if(iplate<1 || iplate>5) {
2017         printf("  iplate<1 or iplate>5, iplate=%i\n",iplate);
2018         return;
2019       }
2020
2021       Float_t nStripes=0;
2022       if(iplate>1) {
2023         for (Int_t i=0;i<iplate-1;i++) {
2024           nStripes += zStrips[i];
2025         }
2026       }
2027
2028       istrip=(Int_t)((zNoise-zLen[iplate-1])/((zLen[iplate]-zLen[iplate-1])/zStrips[iplate-1])); //the strip number in the plate
2029       istrip++;
2030
2031       ipadAlongX = (Int_t)(AliTOFGeometry::NpadX()*gRandom->Rndm())+1;
2032       ipadAlongZ = (Int_t)(AliTOFGeometry::NpadZ()*gRandom->Rndm())+1;
2033       ipad=(Int_t)(ipadAlongZ-1)*AliTOFGeometry::NpadX()+ipadAlongX;    //the pad number
2034
2035       Float_t xStrip=(ipadAlongX-1)*AliTOFGeometry::XPad()+AliTOFGeometry::XPad()*gRandom->Rndm()-0.5*AliTOFGeometry::NpadX()*AliTOFGeometry::XPad();//x-coor.in the strip frame
2036       Float_t zStrip=(ipadAlongZ-1)*AliTOFGeometry::ZPad()+AliTOFGeometry::ZPad()*gRandom->Rndm()-0.5*AliTOFGeometry::NpadZ()*AliTOFGeometry::ZPad();//z-coor.in the strip frame
2037
2038       Int_t nPlace[4]={0,0,0,0};
2039       nPlace[0]=ipad;
2040
2041       Int_t   nActivatedPads=0;
2042       Int_t   nFiredPads=0;
2043       Bool_t  isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE};
2044       Float_t tofAfterSimul[4]={0.,0.,0.,0.};
2045       Float_t qInduced[4]={0.,0.,0.,0.};
2046       Float_t averageTime=0.;
2047       Float_t toffornoise=10.+noiseTof->GetRandom(); // 10 ns offset + parameterization [ns]
2048
2049       BorderEffect(zStrip,xStrip,toffornoise,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect
2050
2051       if(nFiredPads) {
2052         for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
2053           if(isFired[indexOfPad]){// the pad has fired
2054
2055             if(indexOfPad==0) {// the hit belongs to a fired pad
2056               isHitOnFiredPad++;
2057               hitArray[isHitOnFiredPad-1].SetX(0.);
2058               hitArray[isHitOnFiredPad-1].SetY(0.);
2059               hitArray[isHitOnFiredPad-1].SetZ(zNoise);
2060               hitArray[isHitOnFiredPad-1].SetNoise(1);
2061             } // close if(indexOfPad==0)
2062
2063             ipad = nPlace[indexOfPad];
2064
2065             Int_t iMapValue=MapPixels[AliTOFGeometry::NSectors()*(iplate-1)+isector-1][istrip-1][ipad-1];
2066
2067             if(iMapValue==0) {
2068               ipixel++;
2069               if(indexOfPad) nPixEdge++;
2070               MapPixels[AliTOFGeometry::NSectors()*(iplate-1)+isector-1][istrip-1][ipad-1]=ipixel;
2071               pixelArray[ipixel-1].SetGeom(isector,iplate,istrip,ipad);
2072               pixelArray[ipixel-1].AddState(1);
2073               pixelArray[ipixel-1].SetRealTime(tofAfterSimul[indexOfPad]);
2074               pixelArray[ipixel-1].SetHit(isHitOnFiredPad);
2075             } else if( tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) {
2076               pixelArray[iMapValue-1].SetTrack(-1);
2077               pixelArray[iMapValue-1].AddState(1);
2078               pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]);
2079               pixelArray[iMapValue-1].SetHit(isHitOnFiredPad);
2080             }  //end of if(iMapValue==0)
2081
2082           }// close if(isFired[indexOfPad])
2083         }  //end loop on activated pads indexOfPad
2084       } // close if(nFiredPads)
2085     }  //end of NOISE cycle
2086   }
2087
2088   // free used memory
2089   if (fZnoise)
2090     {
2091       delete fZnoise;
2092       fZnoise = 0;
2093     }
2094
2095   if (noiseTof)
2096     {
2097       delete noiseTof;
2098        noiseTof = 0;
2099     }
2100
2101   Int_t nNoiseSignals=0;
2102   Int_t nAll=0;
2103   for(Int_t idummy=1; idummy<ipixel+1; idummy++) {
2104     if(hitArray[pixelArray[idummy-1].GetHit()-1].GetNoise()==1) {
2105       nNoiseSignals++;
2106       if(pixelArray[idummy-1].GetState()>10) nAll++;
2107     }
2108   }
2109
2110   if(fdbg) {
2111     cout << " after adding " << fNoise << " noise hits: " << ipixel << " fired pixels (" << nNoiseSignals << " fired by noise, " << nAll << " noise+track)" << endl;
2112     printf(" %i additional pixels are fired by noise due to edge effect\n",nPixEdge);
2113     cout << " End of adding additional noise hits from outer regions" << endl;
2114   }
2115
2116   Float_t occupancy;
2117   // numberOfPads for AliTOFV4 (Full coverage)
2118   // - to be upgraded checking the used TOF version -
2119   Float_t numberOfPads=AliTOFGeometry::NPadXSector()*AliTOFGeometry::NSectors();
2120   occupancy=100.*ipixel/numberOfPads;   // percentage of fired pads
2121   printf(" Overall TOF occupancy (percentage of fired pads after adding noise) = %f\n",occupancy);
2122   delete [] zLen;
2123   delete [] zStrips;
2124
2125 }
2126
2127
2128 //____________________________________________________________________________
2129 void AliTOFReconstructioner::SetMinDistance(AliTOFRecHit* hitArray, Int_t ilastEntry)
2130 {
2131   //
2132   // Set the distance to the nearest hit for hitArray
2133   // ilastEntry is the index of the last entry of hitArray
2134
2135   // starting the setting for the distance to the nearest TOFhit (cm)
2136   for(Int_t i=0; i<ilastEntry; i++) {
2137
2138     if(hitArray[i].GetFirst()==1 && hitArray[i].GetNoise()==0) { // select the first hit of the track
2139       // hits are not due to noise
2140       Float_t minDistance=10000.,squareDistance; // current values of the (square) distance
2141       Int_t jAtMin=0;                            // index of the hit nearest to the i-th hit
2142       Float_t xhit=hitArray[i].X(); // x coordinate for i-th hit
2143       Float_t yhit=hitArray[i].Y(); // y coordinate for i-th hit
2144       Float_t zhit=hitArray[i].Z(); // z coordinate for i-th hit
2145       //  was    for(Int_t j=0; j<isHitOnFiredPad; j++) {
2146       for(Int_t j=0; j<ilastEntry; j++) {
2147         if(i!=j) {
2148           squareDistance=(hitArray[j].X()-xhit)*(hitArray[j].X()-xhit)+
2149             (hitArray[j].Y()-yhit)*(hitArray[j].Y()-yhit)+
2150             (hitArray[j].Z()-zhit)*(hitArray[j].Z()-zhit);
2151           if(squareDistance<minDistance) {
2152             minDistance=squareDistance;
2153             jAtMin=j;
2154           }
2155         }
2156       }
2157       minDistance=TMath::Sqrt(minDistance);
2158       hitArray[i].SetRmin(minDistance);
2159       if(minDistance==0.) printf(" Rmin=0, i=%i, j=%i, x=%f,y=%f,z=%f\n",i,jAtMin,xhit,yhit,zhit);// it cannot happen
2160     }
2161   }
2162
2163 }
2164
2165 // these lines has to be commented till TPC will provide fPx fPy fPz
2166 // and fL in AliTPChit class
2167 //____________________________________________________________________________
2168 /*
2169 void AliTOFReconstructioner::ReadTPCHits(Int_t ntracks, TTree* treehits, TClonesArray* tpchits, Int_t* iTrackPt, Int_t* iparticle, Float_t* ptTrack, AliTOFTrack* trackArray, Int_t& itrack)
2170 {
2171   //
2172   // Read TPC hits for the current event
2173   //
2174   TParticle *particle=0;
2175   Int_t npions=0;    // number of pions for the current event
2176   Int_t nkaons=0;    // number of kaons for the current event
2177   Int_t nprotons=0;  // number of protons for the current event
2178   Int_t nelectrons=0;// number of electrons for the current event
2179   Int_t nmuons=0;    // number of muons for the current event
2180   Int_t ntotalTPChits=0; // total number of TPC hits for the current event
2181   Int_t idummy=-1;       // dummy var used to count double hit TPC cases
2182   Int_t nTpcDoubleHitsLastRow=0; // number of double TPC hits in the last pad row
2183   Int_t nTpcHitsLastRow=0;       // number of TPC hits in the last pad row
2184   Float_t trdpos[2]={0.,0.};
2185   Float_t pos[3];               // TPC hit position
2186   Float_t mom[3]; // momentum components in the last TPC row
2187   Float_t pt=0., tpclen; // pt: transverse momentum in the last TPC row
2188   Int_t nbytes = 0;
2189   Int_t ipart=0, nhits=0, iprim=0;
2190
2191   itrack=0; // itrack: total number of selected TPC tracks
2192
2193   // speed-up the code
2194   treehits->SetBranchStatus("*",0); // switch off all branches
2195   treehits->SetBranchStatus("TPC*",1); // switch on only TPC
2196
2197   for (Int_t track=0; track<ntracks;track++) {
2198     gAlice->ResetHits();
2199     nbytes += treehits->GetEvent(track);
2200
2201
2202     nhits = tpchits->GetEntriesFast();
2203
2204     for (Int_t hit=0;hit<nhits;hit++) {
2205       ntotalTPChits++;
2206       AliTPChit* tpcHit = (AliTPChit*)tpchits->UncheckedAt(hit);
2207       Int_t row = tpcHit->fPadRow;
2208       ipart    = tpcHit->GetTrack();
2209       if(ipart>=fMaxAllTracks) break;
2210       particle = (TParticle*)gAlice->Particle(ipart);
2211       Int_t pdgCode=particle->GetPdgCode();
2212       // only high momentum tracks
2213       // momentum components at production vertex
2214       Float_t pxvtx = particle->Px();
2215       Float_t pyvtx = particle->Py();
2216       Float_t pzvtx = particle->Pz();
2217       Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx);
2218       if(pvtx>fPBound && row == fLastTPCRow) {
2219         Float_t vx = particle->Vx();
2220         Float_t vy = particle->Vy();
2221         Float_t vr = TMath::Sqrt(vx*vx+vy*vy);
2222         Float_t x = tpcHit->X();
2223         Float_t y = tpcHit->Y();
2224         Float_t z = tpcHit->Z();
2225         pos[0]=x; pos[1]=y; pos[2]=z;
2226
2227         Float_t pxtpc = tpcHit->fPx;
2228         Float_t pytpc = tpcHit->fPy;
2229         Float_t pztpc = tpcHit->fPz;
2230         mom[0]=pxtpc; mom[1]=pytpc; mom[2]=pztpc;
2231         Float_t momtpc = TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc+pztpc*pztpc);
2232
2233         if(x*pxtpc+y*pytpc>0) { // only tracks going out of TPC
2234
2235           Float_t isoutgoing = x*pxtpc+y*pytpc+z*pztpc;
2236           isoutgoing /= (momtpc*TMath::Sqrt(x*x+y*y+z*z));
2237           tpclen = tpcHit->fL;
2238
2239
2240           if(ipart!=idummy) {
2241             if(particle->GetFirstMother() < 0) {
2242               Int_t abspdgCode=TMath::Abs(pdgCode);
2243               switch (abspdgCode) {
2244               case 211:
2245                 npions++;
2246                 break ;
2247               case 321:
2248                 nkaons++;
2249                 break ;
2250               case 2212:
2251                 nprotons++;
2252                 break ;
2253               case 11:
2254                 nelectrons++;
2255                 break ;
2256               case 13:
2257                 nmuons++;
2258                 break ;
2259               }
2260             } // close if(particle->GetFirstMother() < 0)
2261           } // close if(ipart!=idummy)
2262
2263           if(gRandom->Rndm()<fTrackingEfficiency && vr<fRadiusvtxBound && ipart!=idummy) {
2264
2265             itrack++;
2266             if(particle->GetFirstMother() < 0) iprim++;
2267
2268             if(itrack>fMaxTracks) {
2269               cout << "itrack=" << itrack << " > MAXTRACKS=" << fMaxTracks << endl;
2270               return;
2271             } // close if(itrack>fMaxTracks)
2272
2273
2274             iparticle[ipart]=itrack;
2275
2276             trackArray[itrack-1].SetTrack(ipart,pvtx,pdgCode,tpclen,pos,mom,trdpos);
2277
2278             pt=TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc); // pt: transverse momentum at TPC
2279             // Filling iTrackPt[MAXTRACKS] by itrack ordering on Pt
2280             if(itrack==1) {
2281               iTrackPt[itrack-1]=itrack;
2282               ptTrack[itrack-1]=pt;
2283             } else {
2284               for (Int_t i=0; i<itrack-1; i++) {
2285                 if(pt>ptTrack[i]) {
2286                   for(Int_t j=i; j<itrack-1; j++) {
2287                     Int_t k=itrack-1+i-j;
2288                     iTrackPt[k]= iTrackPt[k-1];
2289                     ptTrack[k] = ptTrack[k-1];
2290                   }
2291                   iTrackPt[i]=itrack;
2292                   ptTrack[i]=pt;
2293                   break;
2294                 }
2295                 if(i==itrack-2) {
2296                   iTrackPt[itrack-1]=itrack;
2297                   ptTrack[itrack-1]=pt;
2298                 }
2299               }
2300             }
2301
2302           }  //end of itrack
2303           if(vr>fRadiusvtxBound) nTpcHitsLastRow++;
2304           if(ipart==idummy) nTpcDoubleHitsLastRow++;
2305           idummy=ipart;
2306         }  // close if(x*px+y*py>0)
2307       }  // close if(pvtx>fPBound && row == fLastTPCRow)
2308     }  //end of hits
2309   }  // close loop on tracks
2310
2311
2312   if(fdbg) {
2313     cout << ntotalTPChits << " TPC hits in the last TPC row " << fLastTPCRow << endl;
2314     cout << "   " << nTpcHitsLastRow << " - hits with Rvtx>fRadiusvtxBound=" << fRadiusvtxBound << endl;
2315     cout << "   " << nTpcDoubleHitsLastRow << " double TPC hits" << endl;
2316     cout << itrack    << " - extracted TPC tracks   "     << iprim << " - primary" << endl;
2317     cout << npions    << " primary pions reached TPC"     << endl;
2318     cout << nkaons    << " primary kaons reached TPC"     << endl;
2319     cout << nprotons  << " primary protons reached TPC"   << endl;
2320     cout << nelectrons<< " primary electrons reached TPC" << endl;
2321     cout << nmuons    << " primary muons reached TPC"     << endl;
2322   } // if(fdbg)
2323
2324   Int_t primaryInTPC[6]={0,0,0,0,0,0};
2325   primaryInTPC[0]=npions;
2326   primaryInTPC[1]=nkaons;
2327   primaryInTPC[2]=nprotons;
2328   primaryInTPC[3]=nelectrons;
2329   primaryInTPC[4]=nmuons;
2330   primaryInTPC[5]=npions+nkaons+nprotons+nelectrons+nmuons;
2331
2332   if(fdbg) {
2333     printf("  contents of iTrackPt[MAXTRACKS],PtTrack[MAXTRACKS]\n");
2334     for (Int_t i=0; i<itrack; i++) {
2335       printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]);
2336     }
2337     printf(" Check ordered transverse momentum array\n");
2338     for (Int_t i=itrack-1; i>=0; i--) {
2339       printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]);
2340     }
2341   }// if(fdbg)
2342
2343 }
2344 */
2345 //____________________________________________________________________________
2346 void cylcor(Float_t& x, Float_t& y) {
2347   Float_t rho,phi;
2348
2349   rho=TMath::Sqrt(x*x+y*y);
2350   phi=0.;
2351   if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x);
2352   if(phi<0.) phi=phi+2.*TMath::Pi();
2353   x=rho;
2354   y=phi;
2355
2356 }
2357
2358 //____________________________________________________________________________
2359 void AliTOFReconstructioner::Matching(AliTOFTrack* trackArray, AliTOFRecHit* hitArray, Int_t ***mapPixels, AliTOFPad* pixelArray, Int_t* kTOFhitFirst, Int_t& ipixel, Int_t* iTrackPt, Int_t* iTOFpixel, Int_t ntotTpcTracks)
2360 {
2361   Int_t TestTracks,iTestTrack,itest,wPixel=0,itestc;
2362   Int_t * ntest = new Int_t[fMaxTestTracks];
2363   Int_t * testPixel = new Int_t[fMaxTestTracks];
2364   Float_t wLength=0.,wRho=0.,wZ=0.;
2365   Float_t * testLength = new Float_t[fMaxTestTracks];
2366   Float_t * testRho = new Float_t[fMaxTestTracks];
2367   Float_t * testZ = new Float_t[fMaxTestTracks];
2368   Float_t weight;
2369   Float_t * testWeight = new Float_t[fMaxTestTracks];
2370   Float_t rotationFactor,phi0,coslam,sinlam,helixRadius,xHelixCenter,yHelixCenter,zHelixCenter,helixFactor;
2371   Int_t npixel[5],iMapValue,iwork1,iwork2,iwork3,iwork4,ihit=0;
2372   Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
2373                      1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
2374                      -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
2375                      1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
2376                      1, 1,-1, 0, 1, 1, 2, 0};
2377   Float_t theta0,gpx,gpy,gpz,gp,gpt,gtheta,gx,gy,gz,gr,gxLast,gyLast,gzLast,chargeField;
2378   Float_t sumOfTheta=0.,weightTestTracksOutTof[4];
2379   Float_t s,ds,xRespectToHelixCenter,yRespectToHelixCenter,deltaRadius,fp,xp,yp,grho;
2380   Float_t mass,energy,g;
2381   Int_t itrack=0,itr,particleCharge,istep,iplate=0,iPadAlongX=0;
2382   Int_t itra,t34=0,t32=0,t44=0,t43=0,t42=0;
2383   Int_t wstate=0,m2state=0,wPix;
2384   Int_t idelR=0,idelR1=0,idelR2=0,iRmin=0,iRmin1=0,iRmin2=0;
2385   Float_t massArray[50] = {0.0,0.00051,0.00051,0.0,0.1057,0.1057,0.135,0.1396,0.1396,0.4977,
2386                        0.4936,0.4936,0.9396,0.9383,0.9383,0.4977,0.5488,1.1156,1.1894,1.1926,1.1926,
2387                        1.3149,1.3213,1.6724,0.9396,1.1156,1.1894,1.1926,1.1974,1.3149,
2388                        0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
2389   Float_t delR;
2390   Float_t radius,area,normR,normS,cosAngl;
2391   Int_t iPlateFirst,iTestGmax=0;
2392   Int_t fstate,iPrintM1=0,iPrintM2=0;
2393   Float_t gxExtrap=0.,gyExtrap=0.,gzExtrap=0.;
2394   Float_t avSigZ=0,avSigRPHI=0,avSigP=0,avSigPHI=0,avSigTHETA=0;
2395
2396   Float_t gxW,gyW,gzW;
2397   Float_t length0;
2398   Float_t snr=0;
2399   Int_t indexOfTestTrack;
2400   Float_t zPad,xPad;
2401   Int_t istate=0,imax=0,match,iMaxTestTracksOutTof=0,matchw;
2402   Float_t w,wmax=0.,inverseOfParticleSpeed,w2,smat[9],largestWeightTracksOutTof,sw;
2403   Float_t sumWeightTracksOutTof,sGeomWeigth;
2404   Int_t imatched;
2405   Int_t m10=0,m20=0,m22=0,m23=0;
2406   Int_t PRINT=0;
2407   TParticle *particle;
2408
2409   Float_t time=0.;
2410   itr=ntotTpcTracks;
2411   printf(" itr=%i\n",itr);
2412   for (itra=1; itra<itr+1; itra++) {
2413
2414     Int_t itrack=iTrackPt[itra-1];
2415     if(itrack==0) printf("  iTrackPt[itra-1]=0 for itra=%i\n",itra);
2416     Int_t ipart=trackArray[itrack-1].GetTrack();
2417     Float_t pvtx=trackArray[itrack-1].GetP();
2418     Int_t pdgCode=trackArray[itrack-1].GetPdgCode();
2419     Float_t tpclength=trackArray[itrack-1].GetlTPC();
2420     Float_t x=trackArray[itrack-1].GetRxTPC();
2421     Float_t y=trackArray[itrack-1].GetRyTPC();
2422     Float_t z=trackArray[itrack-1].GetRzTPC();
2423     /* vars used for QA
2424     Float_t RxTPC=x;
2425     Float_t RyTPC=y;
2426     Float_t RzTPC=z;
2427     */
2428     Float_t Wx=x;
2429     Float_t Wy=y;
2430     Float_t Wz=z;
2431     Float_t px=trackArray[itrack-1].GetPxTPC();
2432     Float_t py=trackArray[itrack-1].GetPyTPC();
2433     Float_t pz=trackArray[itrack-1].GetPzTPC();
2434     /* vars used for QA
2435     Float_t pxTPC=px;
2436     Float_t pyTPC=py;
2437     Float_t pzTPC=pz;
2438     */
2439     Float_t p = TMath::Sqrt(px*px+py*py+pz*pz);
2440     /* var used for QA
2441     Float_t pTPC=p;
2442     */
2443     Float_t rho = TMath::Sqrt(x*x+y*y);
2444     Float_t phi=0.;
2445     if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x);
2446     if(phi<0.) phi=phi+2.*TMath::Pi();
2447     /* var used for QA
2448     Float_t phiTPC=phi*kRaddeg;
2449     */
2450     if(fSigmavsp) {
2451       if(p==0) printf(" p=%f in g=0.022/p\n",p);
2452       g=0.022/p;
2453       avSigRPHI += g;      // (cm)
2454       if(rho==0) printf(" rho=%f in phi += g*gRandom->Gaus()/rho\n",rho);
2455       phi += g*gRandom->Gaus()/rho;
2456     } else {
2457       if(rho==0) printf(" rho=%f in phi += (SIGMARPHI*gRandom->Gaus()/rho\n",rho);
2458       phi += (fSigmarphi*gRandom->Gaus()/rho);
2459     }
2460     x=rho*TMath::Cos(phi);
2461     y=rho*TMath::Sin(phi);
2462     /* var used for QA
2463     Float_t zTPC=z;
2464     */
2465     if(fSigmavsp) {
2466       if(p==0) printf(" p=%f in g=0.0275/p\n",p);
2467       g=0.0275/p;
2468       avSigZ += g;      // (cm)
2469       z += g*gRandom->Gaus();
2470     } else {
2471       z += fSigmaZ*gRandom->Gaus();
2472     }
2473
2474     // smearing on TPC momentum
2475
2476     {
2477       Float_t pmom,phi,theta,arg;
2478
2479       pmom=TMath::Sqrt(px*px+py*py+pz*pz);
2480       phi=0.;
2481       if(TMath::Abs(px)>0. || TMath::Abs(py)>0.) phi=TMath::ATan2(py,px);
2482       if(phi<0.) phi=phi+2*TMath::Pi();
2483       arg=1.;
2484       if(pmom>0.) arg=pz/pmom;
2485       theta=0.;
2486       if(TMath::Abs(arg)<=1.) theta=TMath::ACos(arg);
2487
2488       if(fSigmavsp) {
2489         if(pmom<=0) printf(" pmom=%f in g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7\n",pmom);
2490         g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7;
2491         g = 0.01*(g*g*g+1.5)*1.24;
2492         avSigP += g;
2493         pmom *= (1+g*gRandom->Gaus());
2494
2495         if(p<10) {
2496           if(pmom<=0) printf(" pmom=%f in g = 1-TMath::Log(pmom)/TMath::Log(10)\n",pmom);
2497           g = 1-TMath::Log(pmom)/TMath::Log(10);
2498           g = 0.001*(g*g*g+0.3)*0.65;  // (radian)
2499         } else {
2500           g = 0.001*0.3*0.65;
2501         }
2502         avSigPHI += g;
2503         phi += g*gRandom->Gaus();
2504         avSigTHETA += g;
2505         theta += g*gRandom->Gaus();
2506
2507       } else {
2508         pmom *= (1+fSigmap*gRandom->Gaus());
2509         phi += fSigmaPhi*gRandom->Gaus();
2510         theta += fSigmaTheta*gRandom->Gaus();
2511       }
2512       gxW=px;
2513       gyW=py;
2514       gzW=pz;
2515
2516       px=pmom*TMath::Sin(theta)*TMath::Cos(phi);
2517       py=pmom*TMath::Sin(theta)*TMath::Sin(phi);
2518       pz=pmom*TMath::Cos(theta);
2519
2520
2521       if(x*px+y*py<=0) {
2522         x=Wx;
2523         y=Wy;
2524         z=Wz;
2525         px=gxW;
2526         py=gyW;
2527         pz=gzW;
2528       }// if(x*px+y*py<=0)
2529     }
2530
2531     p = TMath::Sqrt(px*px+py*py+pz*pz);
2532
2533     particleCharge=charge[PDGtoGeantCode(pdgCode)-1];
2534     mass=massArray[PDGtoGeantCode(pdgCode)-1];
2535     mass=massArray[8-1];       //we take pion mass for all tracks
2536     //             mass=massArray[14-1];       //here we take proton mass for all tracks
2537     energy=TMath::Sqrt(p*p+mass*mass);
2538     chargeField=particleCharge*fField;
2539
2540     g=fRadLenTPC/( (x*px+y*py)/(rho*p) );
2541
2542     if(g<=0) printf(" error, g<=0: g=%f, itra=%i, x,y,px,py=%f, %f, %f, %f\n",g,itra,x,y,px,py);
2543
2544     theta0=13.6*0.001*TMath::Sqrt(g)*(1.+0.038*TMath::Log(g))*energy/(p*p);
2545
2546
2547     // start Loop on test tracks
2548     sumOfTheta=0.;
2549     for(Int_t i=0;i<4;i++) {
2550       weightTestTracksOutTof[i]=0.;
2551     }
2552
2553     itest=0;
2554     for(Int_t i=0;i<fMaxTestTracks;i++) {
2555       ntest[i]=0;
2556       testPixel[i]=0;
2557       testLength[i]=0.;
2558       testRho[i]=0.;
2559       testZ[i]=0.;
2560       testWeight[i]=0.;
2561     }
2562
2563     iPlateFirst=0;
2564     TestTracks=0;
2565     iTestTrack=0;
2566     iTestGmax=0;
2567
2568     length0=0;
2569
2570     for (indexOfTestTrack=0; indexOfTestTrack<fMaxTestTracks; indexOfTestTrack++) {
2571
2572       iTestTrack++;
2573       gpx=px;
2574       gpy=py;
2575       gpz=pz;
2576       gp=p;
2577       if(indexOfTestTrack) {
2578         gtheta=theta0;
2579         EpMulScatt(gpx,gpy,gpz,gp,gtheta);
2580
2581       } else {
2582         gtheta=0;
2583       }
2584
2585       weight=TMath::Exp(-gtheta*gtheta/(2*theta0*theta0));
2586       sumOfTheta += gtheta;
2587
2588       //    ==========================================================
2589       // Calculate crossing of the track in magnetic field with cylidrical surface
2590       // of radius RTOFINNER
2591       //   chargeField = qB, where q is a charge of a particle in units of e,
2592       //                     B is magnetic field in tesla
2593       //   see 3.3.1.1. in the book "Data analysis techniques for
2594       //   high-energy physics experiments", edited by M.Regler
2595       //   in Russian: "Metody analiza dannykh v fizicheskom eksperimente"
2596       //   Moskva, "Mir", 1993. ctr.306
2597
2598       // Initial constants
2599       rotationFactor=1.;
2600       if(chargeField<0.) rotationFactor=-1.;
2601       rotationFactor=-rotationFactor;
2602       gpt=gpx;
2603       phi0=gpy;
2604       cylcor(gpt,phi0);
2605       phi0 -= rotationFactor*TMath::Pi()*0.5;
2606       //               phi0 -= h*PID2;
2607       coslam=gpt/gp;
2608       sinlam=gpz/gp;
2609       //      helixRadius=100.*gpt/TMath::Abs(0.299792458*chargeField);
2610       helixRadius=100.*gpt/TMath::Abs(AliTOFGeometry::SpeedOfLight()*chargeField);
2611       xHelixCenter=x-helixRadius*TMath::Cos(phi0);
2612       yHelixCenter=y-helixRadius*TMath::Sin(phi0);
2613       zHelixCenter=z;
2614       helixFactor=rotationFactor*coslam/helixRadius;
2615
2616       //   Solves the equation f(s)=r(s)-RTOFINNER=0 by the Newton's method:
2617       //   snew=s-f/f'
2618       istep=0;
2619       s=AliTOFGeometry::Rmin()-TMath::Sqrt(x*x+y*y);;
2620       do {
2621         istep++;
2622         xRespectToHelixCenter=helixRadius*TMath::Cos(phi0+s*helixFactor);
2623         yRespectToHelixCenter=helixRadius*TMath::Sin(phi0+s*helixFactor);
2624         gx=xHelixCenter+xRespectToHelixCenter;
2625         gy=yHelixCenter+yRespectToHelixCenter;
2626         gr=TMath::Sqrt(gx*gx+gy*gy);
2627         deltaRadius=gr-AliTOFGeometry::Rmin();
2628         xp=-helixFactor*yRespectToHelixCenter;
2629         yp= helixFactor*xRespectToHelixCenter;
2630         fp=(gx*xp+gy*yp)/gr;
2631         ds=deltaRadius/fp;
2632         s -= ds;
2633         if(istep==20) {
2634           istep=0;
2635           break;
2636         }
2637       } while (TMath::Abs(ds)>0.01);
2638
2639
2640       if(istep==0) goto end;
2641
2642       //   Steps along the circle till a pad
2643       wPixel=0;
2644       wLength=0.;
2645       iplate=0;
2646       iPadAlongX=0;
2647       grho=0.;
2648       ds=fStep;
2649       gxLast=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor);
2650       gyLast=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor);
2651       gzLast=zHelixCenter+s*sinlam;
2652
2653
2654       do {
2655         istep++;
2656         s += ds;
2657         gx=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor);
2658         gy=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor);
2659         gz=zHelixCenter+s*sinlam;
2660         rho=TMath::Sqrt(gx*gx+gy*gy);
2661
2662         IsInsideThePad(gMC,gx,gy,gz,npixel,zPad,xPad);
2663
2664         iplate += npixel[1];
2665         iPadAlongX += npixel[4];
2666
2667         if(indexOfTestTrack==0 && iplate && iPlateFirst==0) {
2668           iPlateFirst=1;
2669           length0=s;
2670
2671           radius=s*3*theta0;
2672           area=TMath::Pi()*radius*radius;
2673           normR=TMath::Sqrt(gx*gx+gy*gy);
2674           normS=TMath::Sqrt((gx-gxLast)*(gx-gxLast)+
2675                      (gy-gyLast)*(gy-gyLast)+
2676                      (gz-gzLast)*(gz-gzLast));
2677
2678           cosAngl=(gx*(gx-gxLast)+gy*(gy-gyLast))/(normR*normS);
2679           if(cosAngl<0) printf(" cosAngl<0: gx=%f,gy=%f,  gxLast=%f,gyLast=%f,gzLast=%f\n",gx,gy,gxLast,gyLast,gzLast);
2680
2681           area /= cosAngl;
2682           TestTracks=(Int_t) (2*area/(AliTOFGeometry::XPad() * AliTOFGeometry::ZPad()));
2683
2684           if(TestTracks<12) TestTracks=12;
2685
2686           // Angles of entering into the TOF plate
2687
2688           Int_t iZ=0;
2689           if(TMath::Abs(gz)>300) {
2690             iZ=4;
2691           } else if(TMath::Abs(gz)>200) {
2692             iZ=3;
2693           } else if(TMath::Abs(gz)>100) {
2694             iZ=2;
2695           } else if(TMath::Abs(gz)>0) {
2696             iZ=1;
2697           }
2698
2699
2700         } // end of if(indexOfTestTrack==0 && iplate && iPlateFirst==0)
2701
2702
2703         if(npixel[4]>0) {
2704
2705           iwork1=npixel[0];
2706           iwork2=npixel[1];
2707           iwork3=npixel[2];
2708           //               iwork4=npixel[3];
2709           iwork4=(npixel[3]-1)*AliTOFGeometry::NpadX()+npixel[4];
2710
2711           Int_t ifirstindex=AliTOFGeometry::NSectors()*(npixel[1]-1)+npixel[0];
2712           iMapValue=mapPixels[ifirstindex-1][iwork3-1][iwork4-1];
2713           if(iMapValue==0) {
2714             ipixel++;
2715             if(ipixel>fMaxPixels) {
2716               cout << "ipixel=" << ipixel << " > MAXPIXELS=" << fMaxPixels << endl;
2717               break;
2718             }
2719             mapPixels[ifirstindex-1][iwork3-1][iwork4-1]=ipixel;
2720             pixelArray[ipixel-1].SetGeom(iwork1,iwork2,iwork3,iwork4);
2721             iMapValue=ipixel;
2722           }
2723
2724           wPixel=iMapValue;
2725           wLength=tpclength+s;
2726           wRho=rho;
2727           wZ=gz;
2728
2729           ihit=kTOFhitFirst[ipart];
2730
2731           if(ihit) {
2732             if(indexOfTestTrack==0) {
2733               {
2734                 idelR++;
2735                 delR=TMath::Sqrt((gx-hitArray[ihit-1].X())*(gx-hitArray[ihit-1].X())+
2736                           (gy-hitArray[ihit-1].Y())*(gy-hitArray[ihit-1].Y())+
2737                           (gz-hitArray[ihit-1].Z())*(gz-hitArray[ihit-1].Z()));
2738
2739               }
2740
2741               if(delR>hitArray[ihit-1].GetRmin()) iRmin++;
2742               gxExtrap=gx;
2743               gyExtrap=gy;
2744               gzExtrap=gz;
2745             } else {
2746               delR=TMath::Sqrt((gx-gxExtrap)*(gx-gxExtrap)+
2747                         (gy-gyExtrap)*(gy-gyExtrap)+
2748                         (gz-gzExtrap)*(gz-gzExtrap));
2749             }
2750           }  //end of if(ihit)
2751
2752           break;
2753
2754         }  //end of npixel[4]
2755
2756         if(rho<grho) {
2757           istep=0;
2758           break;
2759         }
2760         grho=rho;
2761
2762         gxLast=gx;
2763         gyLast=gy;
2764         gzLast=gz;
2765
2766       } while(rho<AliTOFGeometry::Rmax()); //end of do
2767
2768
2769       if(istep>0) {
2770         if(iplate) {
2771           if(iPadAlongX==0) {
2772             istep=-3;            // holes in TOF
2773           }
2774         } else {
2775           if(TMath::Abs(gz)<AliTOFGeometry::MaxhZtof()) {
2776             //                   if(TMath::Abs(gz)<MAXZTOF2) {
2777             istep=-2;            // PHOS and RICH holes or holes in between TOF plates
2778           } else {
2779             istep=-1;            // out of TOF on z-size
2780           }
2781         }
2782       }
2783
2784       if(iPadAlongX>0) {
2785         if(itest==0) {
2786           itest=1;
2787           ntest[itest-1]=1;
2788           testPixel[itest-1]=wPixel;
2789           testLength[itest-1]=wLength;
2790           testRho[itest-1]=wRho;
2791           testZ[itest-1]=wZ;
2792           testWeight[itest-1]=weight;
2793         } else {
2794           Int_t k=0;
2795           for(Int_t i=0;i<itest;i++) {
2796             k=0;
2797             if(testPixel[i]==wPixel) {
2798               k=1;
2799               ntest[i]++;
2800               testLength[i] += wLength;
2801               testRho[i] += wRho;
2802               testZ[i] += wZ;
2803               testWeight[i] += weight;
2804               break;
2805             }
2806           }  //end for i
2807           if(k==0) {
2808             itest++;
2809             ntest[itest-1]=1;
2810             testPixel[itest-1]=wPixel;
2811             testLength[itest-1]=wLength;
2812             testRho[itest-1]=wRho;
2813             testZ[itest-1]=wZ;
2814             testWeight[itest-1]=weight;
2815           }
2816         }
2817       }
2818
2819     end: ;
2820       //   Statistics
2821       if(fMatchingStyle==1) {
2822         if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] ++;
2823       } else {
2824         if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] += weight;
2825       }
2826
2827       if(fMatchingStyle==2) {
2828         if(indexOfTestTrack==0 && istep==0) break;
2829         if(indexOfTestTrack+1==TestTracks) break;
2830       }
2831
2832     }  //end of indexOfTestTrack
2833
2834     snr += (Float_t) (indexOfTestTrack+1);
2835
2836     //   Search for the "hole" with the largest weigth
2837     largestWeightTracksOutTof=0.;
2838     sumWeightTracksOutTof=0.;
2839     for(Int_t i=0;i<4;i++) {
2840       w=weightTestTracksOutTof[i];
2841       sumWeightTracksOutTof += w;
2842       if(w>largestWeightTracksOutTof) {
2843         largestWeightTracksOutTof=w;
2844         iMaxTestTracksOutTof=i;
2845       }
2846     }
2847
2848     itestc=itest;
2849     if(itest>0) {
2850       for(Int_t i=0;i<itest;i++) {
2851         testLength[i] /= ntest[i];
2852         testRho[i] /= ntest[i];
2853         testZ[i] /= ntest[i];
2854       }
2855       //   Search for the pixel with the largest weigth
2856       wmax=0.;
2857       wstate=0;
2858       sw=0;
2859       sGeomWeigth=0;
2860       for(Int_t i=0;i<itest;i++) {
2861         istate=pixelArray[testPixel[i]-1].GetState();
2862         fstate=0;
2863         if(istate>0) {
2864           fstate=1;
2865           wstate++;
2866         }
2867         if(fMatchingStyle==1) {
2868           sGeomWeigth += ntest[i];
2869           w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*ntest[i];
2870           if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1;
2871         } else {
2872           sGeomWeigth += testWeight[i];
2873           w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*testWeight[i];
2874           if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1;
2875         }
2876
2877         // weighting according to the Pulse Height (we use the square of weight)
2878         // if (fChargeFactorForMatching) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge());
2879         if (fChargeFactorForMatching && fstate==1) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge());
2880
2881         if(w>wmax) {
2882           wmax=w;
2883           imax=i;
2884         }
2885         sw += w;
2886       }
2887       wPixel=testPixel[imax];
2888       wLength=testLength[imax];
2889       istate=pixelArray[wPixel-1].GetState();
2890
2891       //Choose the TOF dead space
2892       //               if(istate==0 && largestWeightTracksOutTof>wmax) {
2893       //               if(istate==0 && largestWeightTracksOutTof>=sw) {
2894       if(istate==0 && sumWeightTracksOutTof>sGeomWeigth) {
2895         itestc=itest;
2896         itest=0;
2897       }
2898     }
2899
2900     if(itest>0) {
2901
2902       //   Set for MyTrack: Pixel
2903       trackArray[itrack-1].SetPixel(wPixel);
2904
2905       istate=pixelArray[wPixel-1].GetState();
2906
2907       if(istate) {
2908
2909         //   Set for MyTrack: Pixel, Length, TOF, MassTOF
2910         //fp
2911         //time=pixelArray[wPixel-1].GetTime();
2912         time=pixelArray[wPixel-1].GetRealTime();
2913         trackArray[itrack-1].SetLength(wLength);
2914         trackArray[itrack-1].SetTof(time);
2915
2916         inverseOfParticleSpeed=time/wLength;
2917         //w=900.*inverseOfParticleSpeed*inverseOfParticleSpeed-1.;
2918         w=(100.*AliTOFGeometry::SpeedOfLight())*(100.*AliTOFGeometry::SpeedOfLight())*inverseOfParticleSpeed*inverseOfParticleSpeed-1.;
2919         w2=pvtx*pvtx;
2920         Float_t squareMass=w2*w;
2921         mass=TMath::Sqrt(TMath::Abs(squareMass));
2922         if(w<0.) mass=-mass;
2923
2924         trackArray[itrack-1].SetMassTOF(mass);
2925
2926         //   Set for MyTrack: Matching
2927         match=4;
2928         //                 if(ipart==pixelArray[wPixel-1].GetTrack()) match=3;
2929         if( (ipart==pixelArray[wPixel-1].GetTrack()) && hitArray[pixelArray[wPixel-1].GetHit()-1].GetNoise()==0)match=3;
2930         imatched=pixelArray[wPixel-1].GetTrackMatched();
2931         //   Set for TOFPixel the number of matched track
2932         pixelArray[wPixel-1].SetTrackMatched(itrack);
2933
2934         if(imatched>0) {
2935           matchw=trackArray[imatched-1].GetMatching();
2936           if(match==3 && matchw==4) t34++;
2937           if(match==3 && matchw==2) t32++;
2938           if(match==4 && matchw==4) t44++;
2939           if(match==4 && matchw==3) t43++;
2940           if(match==4 && matchw==2) t42++;
2941           if(iTOFpixel[ipart]==0 || iTOFpixel[trackArray[imatched-1].GetTrack()]==0) {
2942             m20++;
2943           } else if(iTOFpixel[ipart]==iTOFpixel[trackArray[imatched-1].GetTrack()]) {
2944             m22++;
2945           } else {
2946             m23++;
2947             wPix=iTOFpixel[ipart];
2948             if(PRINT && iPrintM1==10 && iPrintM2<10) {
2949               if(iPrintM2==0) {
2950                 printf("*** test print for tracks matched with the pixel for with we had matched track\n");
2951               }
2952               iPrintM2++;
2953               printf(" m=2: ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n",
2954                      ipart,pdgCode,p,theta0,wPix,
2955                      pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel());
2956               printf("      mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n",
2957                      match,wPixel,
2958                      pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(),
2959                      itest,imax,wmax,testZ[imax],wstate);
2960               Int_t fstat,istat;
2961               for(Int_t i=0;i<itest;i++) {
2962                 wPix=testPixel[i];
2963                 istat=pixelArray[wPix-1].GetState();
2964                 fstat=0;
2965                 if(istat>0) fstat=1;
2966                 w=(fpadefficiency*fstat+(1.-fpadefficiency)*(1-fstat))*ntest[i];
2967                 if(istat>0)
2968                   printf("                     %i: %i Pixel(LP=%i,SP=%i,P=%i), istat=%i, ntest=%i, w=%f\n",i+1,
2969                          wPix,pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel(),
2970                          istat,ntest[i],w);
2971               }
2972               printf("      mat=%i, %i Pixel \n",matchw,trackArray[imatched-1].GetPad());
2973             }
2974           }
2975           if(wstate>1) m2state++;
2976           smat[matchw+4]--;
2977           match=2;
2978           trackArray[imatched-1].SetMatching(match);
2979           smat[match+4]++;
2980
2981         }  // if(imatched>0)
2982
2983       } else {  //else if(istate)
2984
2985         match=1;
2986         if(iTOFpixel[ipart]==0) m10++;
2987         if(PRINT && iPrintM1<10) {
2988           Int_t wPix;
2989           wPix=iTOFpixel[ipart];
2990           if(wPix) {
2991             if(iPrintM1==0) {
2992               printf("*** test print for tracks fired a pixel but matched with non-fired pixel\n");
2993             }
2994             iPrintM1++;
2995             printf(" m=1: itra=%i,ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n",
2996                    itra,ipart,pdgCode,p,theta0,wPix,
2997                    pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel());
2998             printf("      mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n",
2999                    match,wPixel,
3000                    pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(),
3001                    itest,imax,wmax,testZ[imax],wstate);
3002
3003           }
3004         }  //end if(PRINT && iPrintM1<10)
3005
3006       }  //end if(istate)
3007
3008     } else {
3009       match=-1-iMaxTestTracksOutTof;
3010
3011     }  //end itest
3012
3013     trackArray[itrack-1].SetMatching(match);
3014     //             if(iTestGmax==1) hMTT->Fill(match);
3015     smat[match+4]++;
3016
3017     sumOfTheta /= iTestTrack;
3018
3019     itest=itestc;
3020
3021     //Test
3022     if(PRINT) {
3023       if(iTOFpixel[ipart] && match!=3) {
3024         particle = (TParticle*)gAlice->GetMCApp()->Particle(ipart);  //for V3.05
3025
3026         printf("          ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), fired by %i track\n",iTOFpixel[ipart],
3027                         pixelArray[iTOFpixel[ipart]-1].GetSector(),pixelArray[iTOFpixel[ipart]-1].GetPlate(),
3028                         pixelArray[iTOFpixel[ipart]-1].GetStrip(),pixelArray[iTOFpixel[ipart]-1].GetPixel(),
3029                         pixelArray[iTOFpixel[ipart]-1].GetTrack());
3030         printf("     indexOfTestTrack=%i itest=%i weightTestTracksOutTof[4]=%f weightTestTracksOutTof[2]=%f weightTestTracksOutTof[1]=%f weightTestTracksOutTof[0]=%f\n",
3031                         indexOfTestTrack,itest,weightTestTracksOutTof[3],weightTestTracksOutTof[2],
3032                         weightTestTracksOutTof[1],weightTestTracksOutTof[0]);
3033         if(itest) {
3034
3035           printf("     take ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), (fired by %i track), match=%i\n",
3036                           wPixel,pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetStrip(),
3037                           pixelArray[wPixel-1].GetPixel(),pixelArray[wPixel-1].GetTrack(),match);
3038         }
3039       }
3040     }
3041     if(PRINT && itra<10 ) {
3042
3043       if(itest) {
3044         cout << "      number of pixels with test tracks=" << itest << endl;
3045         for(Int_t i=0;i<itest;i++) {
3046           cout << "      " << i+1 << "  tr.=" << ntest[i] << "  w=" << testWeight[i] << "  pix.= " << testPixel[i] << " (" <<
3047             pixelArray[testPixel[i]-1].GetSector() << " " << " " << pixelArray[testPixel[i]-1].GetPlate() << " " <<
3048             pixelArray[testPixel[i]-1].GetPixel() << " )" << "  l= " << testLength[i] << " sig=" <<
3049             theta0*(testLength[i]-tpclength) << "  rho= " << testRho[i] << "  z= " << testZ[i] << endl;
3050         }
3051         cout << "      pixel=" << wPixel << "  state=" << istate << "  l=" << wLength << "  TOF=" << time << "  m=" << mass << "  match=" << match <<  endl;
3052         if(istate>0) cout << "      fired by track " << pixelArray[wPixel-1].GetTrack() << endl;
3053       }
3054     }
3055   }  //end of track
3056
3057
3058   if(itr) {
3059     printf(" %f probe tracks per 1 real track\n",snr/itr);
3060     itrack=itr;
3061   }
3062
3063
3064   cout << ipixel << " - total number of TOF pixels after matching" << endl;
3065   w=iRmin;
3066   if(idelR!=0) {
3067     w /= idelR;
3068     printf(" %i tracks with delR, %f of them have delR>Rmin \n",idelR,w);
3069   }
3070   w=iRmin1;
3071   if(idelR1!=0) {
3072     w /= idelR1;
3073     printf(" %i tracks with delR1 (|z|<175), %f of them have delR>Rmin \n",idelR1,w);
3074   }
3075   w=iRmin2;
3076   if(idelR2!=0) {
3077     w /= idelR2;
3078     printf(" %i tracks with delR2 (|z|>175), %f of them have delR>Rmin \n",idelR2,w);
3079   }
3080
3081   cout << " ********************  End of matching **********" << endl;
3082   delete [] ntest;
3083   delete [] testPixel;
3084   delete [] testLength;
3085   delete [] testRho;
3086   delete [] testZ;
3087   delete [] testWeight;
3088 }
3089
3090 //____________________________________________________________________________
3091 void AliTOFReconstructioner::FillNtuple(Int_t ntracks, AliTOFTrack* trackArray, AliTOFRecHit* hitArray, AliTOFPad* pixelArray, Int_t* iTOFpixel, Int_t* iparticle, Float_t* toftime, Int_t& ipixelLastEntry, Int_t itrack){
3092
3093   // itrack : total number of TPC selected tracks
3094   // for the caller is ntotTPCtracks
3095
3096   cout << " ********************  Start of searching non-matched fired pixels **********" << endl;
3097   const Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
3098                            1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
3099                            -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
3100                            1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
3101                            1, 1,-1, 0, 1, 1, 2, 0};
3102
3103   Int_t macthm1=0;
3104   Int_t macthm2=0;
3105   Int_t macthm3=0;
3106   Int_t macthm4=0;
3107   Int_t macth0=0;
3108   Int_t macth1=0;
3109   Int_t macth2=0;
3110   Int_t macth3=0;
3111   Int_t macth4=0;
3112
3113
3114   Float_t smat[9],smat0[9],smat1[9];
3115   for(Int_t i=0;i<9;i++) {
3116     smat[i]=0.;
3117     smat0[i]=0.;
3118     smat1[i]=0.;
3119   }
3120
3121   Int_t nFiredPixelsNotMatchedWithTracks=0;
3122   Int_t istate;
3123   for (Int_t i=0; i<ipixelLastEntry; i++) {
3124     istate=pixelArray[i].GetState();
3125     if(istate==0) break;
3126     if(pixelArray[i].GetTrackMatched()==-1) nFiredPixelsNotMatchedWithTracks++;
3127   }
3128   printf("  %i fired pixels have not matched tracks\n",nFiredPixelsNotMatchedWithTracks);
3129   cout << " ********************  End of searching non-matched fired pixels **********" << endl;
3130
3131   Int_t nTPCHitMissing=0;
3132   for(Int_t i=0; i<ipixelLastEntry; i++) {
3133     if(pixelArray[i].GetHit()>0) {
3134       if(hitArray[pixelArray[i].GetHit()-1].GetNoise()==0) {
3135         if(iparticle[pixelArray[i].GetTrack()]==0) nTPCHitMissing++;
3136       }
3137     }
3138   }
3139   printf("  %i pixels fired by track hit without a hit on the last layer of TPC\n",nTPCHitMissing);
3140
3141
3142   Int_t icharge=0;   // total number of charged particles
3143   Int_t iprim=0;     // number of primaries
3144   Int_t ipions=0;    // number of primary pions
3145   Int_t ikaons=0;    // number of primary kaons
3146   Int_t iprotons=0;  // number of primary protons
3147   Int_t ielectrons=0;// number of primary electrons
3148   Int_t imuons=0;    // number of primary muons
3149   Float_t particleTypeArray[6][5][2];
3150
3151   for (Int_t index1=0;index1<6;index1++) {
3152     for (Int_t index2=0;index2<5;index2++) {
3153       for (Int_t index3=0;index3<2;index3++) {
3154         particleTypeArray[index1][index2][index3]=0.;
3155       }
3156     }
3157   }
3158
3159   Int_t nTOFhitsWithNoTPCTracks=0; // to be moved later when used
3160
3161   /*
3162   TObjArray *Particles = gAlice->Particles();
3163   Int_t numberOfParticles=Particles->GetEntries();
3164   cout << "numberOfParticles " << numberOfParticles << endl;
3165   // fpdbg
3166   if(numberOfParticles>fMaxAllTracks) numberOfParticles=fMaxAllTracks;
3167   */
3168
3169   for (Int_t i=0; i<ntracks; i++) { // starting loop on all primaries charged particles for current event)
3170
3171     /*
3172     cout << "particle " << i << endl;
3173     cout << "total " << numberOfParticles << endl;
3174     */
3175     TParticle *part = (TParticle *) gAlice->GetMCApp()->Particle(i);
3176     if(charge[PDGtoGeantCode(part->GetPdgCode())-1]) {
3177       icharge++;
3178       /*
3179       cout << "charged particles " << icharge << endl;
3180       */
3181       Int_t particleType=0;
3182       Int_t absPdgCode = TMath::Abs(part->GetPdgCode());
3183       switch (absPdgCode) {
3184       case 211:
3185         particleType=3;
3186         break ;
3187       case 321:
3188         particleType=2;
3189         break ;
3190       case 2212:
3191         particleType=1;
3192         break ;
3193       case 11:
3194         particleType=4;
3195         break ;
3196       case 13:
3197         particleType=5;
3198         break ;
3199       }
3200
3201       if(part->GetFirstMother() < 0) {
3202         iprim++;
3203         switch (particleType) {
3204         case 1:
3205           iprotons++;
3206           break ;
3207         case 2:
3208           ikaons++;
3209           break ;
3210         case 3:
3211           ipions++;
3212           break ;
3213         case 4:
3214           ielectrons++;
3215           break ;
3216         case 5:
3217           imuons++;
3218           break ;
3219         }
3220       }
3221
3222       Int_t match=0;
3223       Float_t wLength=-1.;
3224       Float_t time=-1.;
3225       Float_t mass=-1.;
3226
3227       Int_t itr=iparticle[i]; // get the track number for the current charged particle
3228
3229       if(iTOFpixel[i]>0 && itr==0) nTOFhitsWithNoTPCTracks++;
3230
3231       if(itr) {
3232         match=trackArray[itr-1].GetMatching();
3233         //cout << "match " << match << endl;
3234         wLength=trackArray[itr-1].GetLength();
3235         //cout << "wLength " << wLength << endl;
3236         time=trackArray[itr-1].GetTof();
3237         mass=trackArray[itr-1].GetMassTOF();
3238         //cout << "mext " << mass << endl;
3239         //        if(PRINT && (i>789 && i<800) ) cout << i << " track:  l=" << wLength << "  TOF=" << time << "  m=" << mass << "  match=" << match <<  endl;
3240         if(iTOFpixel[i]==0) {
3241           smat0[match+4]++;