TOF/AliTOFReconstructioner.cxx

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