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