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