TOF/AliTOFSDigitizer.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 // This is a TTask that constructs SDigits out of Hits
  20 // A Summable Digits is the "sum" of all hits in a pad
  21 // Detector response has been simulated via the method
  22 // SimulateDetectorResponse
  23 //
  24 //-- Authors: F. Pierella, A. De Caro
  25 // Use case: see AliTOFhits2sdigits.C macro in the CVS
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28 #include <Riostream.h>
  29 #include <stdlib.h>
  30
  31 #include <TBenchmark.h>
  32 #include <TF1.h>
  33 #include <TFile.h>
  34 #include <TFolder.h>
  35 #include <TH1.h>
  36 #include <TParticle.h>
  37 #include <TROOT.h>
  38 #include <TSystem.h>
  39 #include <TTask.h>
  40 #include <TTree.h>
  41
  42 #include "AliLog.h"
  43 #include "AliDetector.h"
  44 #include "AliLoader.h"
  45 #include "AliMC.h"
  46 #include "AliRun.h"
  47 #include "AliRunLoader.h"
  48
  49 #include "AliTOF.h"
  50 #include "AliTOFGeometry.h"
  51 #include "AliTOFHitMap.h"
  52 #include "AliTOFSDigit.h"
  53 #include "AliTOFSDigitizer.h"
  54 #include "AliTOFhit.h"
  55 #include "AliTOFhitT0.h"
  56
  57 ClassImp(AliTOFSDigitizer)
  58
  59 //____________________________________________________________________________
  60   AliTOFSDigitizer::AliTOFSDigitizer():TTask("TOFSDigitizer","")
  61 {
  62   // ctor
  63
  64   fRunLoader      = 0;
  65   fTOFLoader      = 0;
  66
  67   fEvent1         = 0;
  68   fEvent2         = 0;
  69   ftail           = 0;
  70   fSelectedSector = -1;
  71   fSelectedPlate  = -1;
  72
  73   fTOFGeometry = new AliTOFGeometry();
  74
  75 }
  76
  77 //------------------------------------------------------------------------
  78 AliTOFSDigitizer::AliTOFSDigitizer(const AliTOFSDigitizer &source)
  79   :TTask(source)
  80 {
  81   // copy constructor
  82   this->fTOFGeometry=source.fTOFGeometry;
  83
  84 }
  85
  86 //____________________________________________________________________________
  87 AliTOFSDigitizer& AliTOFSDigitizer::operator=(const AliTOFSDigitizer &source)
  88 {
  89   // ass. op.
  90   this->fTOFGeometry=source.fTOFGeometry;
  91   return *this;
  92
  93 }
  94
  95 //____________________________________________________________________________
  96 AliTOFSDigitizer::AliTOFSDigitizer(const char* HeaderFile, Int_t evNumber1, Int_t nEvents):TTask("TOFSDigitizer","")
  97 {
  98   //ctor, reading from input file
  99   ftail    = 0;
 100   fSelectedSector=-1; // by default we sdigitize all sectors
 101   fSelectedPlate =-1; // by default we sdigitize all plates in all sectors
 102
 103   fHeadersFile = HeaderFile ; // input filename (with hits)
 104   TFile * file = (TFile*) gROOT->GetFile(fHeadersFile.Data());
 105
 106   //File was not opened yet open file and get alirun object
 107   if (file == 0) {
 108     file   = TFile::Open(fHeadersFile.Data(),"update") ;
 109     gAlice = (AliRun *) file->Get("gAlice") ;
 110   }
 111
 112   // add Task to //root/Tasks folder
 113   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 114   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
 115   if (!fRunLoader)
 116     fRunLoader = AliRunLoader::Open(HeaderFile);//open session and mount on default event folder
 117   if (fRunLoader == 0x0)
 118     {
 119       AliFatal("Event is not loaded. Exiting");
 120       return;
 121     }
 122
 123   fRunLoader->CdGAFile();
 124   TDirectory *savedir=gDirectory;
 125   TFile *in=(TFile*)gFile;
 126
 127   if (!in->IsOpen()) {
 128     AliWarning("Geometry file is not open default TOF geometry will be used");
 129     fTOFGeometry = new AliTOFGeometry();
 130   }
 131   else {
 132     in->cd();
 133     fTOFGeometry = (AliTOFGeometry*)in->Get("TOFgeometry");
 134   }
 135
 136   savedir->cd();
 137
 138   if (fRunLoader->TreeE() == 0x0) fRunLoader->LoadHeader();
 139
 140   if (evNumber1>=0) fEvent1 = evNumber1;
 141   else fEvent1=0;
 142
 143   if (nEvents==0) fEvent2 = (Int_t)(fRunLoader->GetNumberOfEvents());
 144   else if (nEvents>0) fEvent2 = evNumber1+nEvents;
 145   else fEvent2 = 1;
 146
 147   if (!(fEvent2>fEvent1)) {
 148     AliError(Form("fEvent2 = %d <= fEvent1 = %d", fEvent2, fEvent1));
 149     fEvent1 = 0;
 150     fEvent2 = 1;
 151     AliError(Form("Correction: fEvent2 = %d <= fEvent1 = %d", fEvent2, fEvent1));
 152   }
 153
 154   // init parameters for sdigitization
 155   InitParameters();
 156
 157   fTOFLoader = fRunLoader->GetLoader("TOFLoader");
 158   if (fTOFLoader == 0x0)
 159     {
 160       AliFatal("Can not find TOF loader in event. Exiting.");
 161       return;
 162     }
 163   fTOFLoader->PostSDigitizer(this);
 164 }
 165
 166 //____________________________________________________________________________
 167 AliTOFSDigitizer::~AliTOFSDigitizer()
 168 {
 169   // dtor
 170   fTOFLoader->CleanSDigitizer();
 171
 172   delete fTOFGeometry;
 173
 174 }
 175
 176 //____________________________________________________________________________
 177 void AliTOFSDigitizer::InitParameters()
 178 {
 179   // set parameters for detector simulation
 180
 181   fTimeResolution = 0.080; //0.120; OLD
 182   fpadefficiency  = 0.99 ;
 183   fEdgeEffect     = 2   ;
 184   fEdgeTails      = 0   ;
 185   fHparameter     = 0.4 ;
 186   fH2parameter    = 0.15;
 187   fKparameter     = 0.5 ;
 188   fK2parameter    = 0.35;
 189   fEffCenter      = fpadefficiency;
 190   fEffBoundary    = 0.65;
 191   fEff2Boundary   = 0.90;
 192   fEff3Boundary   = 0.08;
 193   fAddTRes        = 68. ; // \sqrt{2x20^2 + 15^2 + 2x10^2 + 30^2 + 50^2} (p-p)
 194   //fAddTRes      = 48. ; // \sqrt{2x20^2 + 15^2 + 2x10^2 + 30^2 + 15^2} (Pb-Pb)
 195   // 30^2+20^2+40^2+50^2+50^2+50^2 = 10400 ps^2 (very old value)
 196   fResCenter      = 35. ; //50. ; // OLD
 197   fResBoundary    = 70. ;
 198   fResSlope       = 37. ; //40. ; // OLD
 199   fTimeWalkCenter = 0.  ;
 200   fTimeWalkBoundary=0.  ;
 201   fTimeWalkSlope  = 0.  ;
 202   fTimeDelayFlag  = 1   ;
 203   fPulseHeightSlope=2.0 ;
 204   fTimeDelaySlope =0.060;
 205   // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.);
 206   fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter);
 207   fChargeSmearing=0.0   ;
 208   fLogChargeSmearing=0.13;
 209   fTimeSmearing   =0.022;
 210   fAverageTimeFlag=0    ;
 211
 212   fAdcBin   = 0.25;    // 1 ADC bin = 0.25 pC (or 0.03 pC)
 213   fAdcMean  = 50.;     // ADC distribution mpv value for Landau (in bins)
 214                        // it corresponds to a mean value of ~100 bins
 215   fAdcRms   = 25.;     // ADC distribution rms value (in bins)
 216                        // it corresponds to distribution rms ~50 bins
 217 }
 218
 219 //__________________________________________________________________
 220 Double_t TimeWithTail(Double_t* x, Double_t* par)
 221 {
 222   // sigma - par[0], alpha - par[1], part - par[2]
 223   //  at x<part*sigma - gauss
 224   //  at x>part*sigma - TMath::Exp(-x/alpha)
 225   Float_t xx =x[0];
 226   Double_t f;
 227   if(xx<par[0]*par[2]) {
 228     f = TMath::Exp(-xx*xx/(2*par[0]*par[0]));
 229   } else {
 230     f = TMath::Exp(-(xx-par[0]*par[2])/par[1]-0.5*par[2]*par[2]);
 231   }
 232   return f;
 233 }
 234
 235 //____________________________________________________________________________
 236 void AliTOFSDigitizer::Exec(Option_t *verboseOption) {
 237   //execute TOF sdigitization
 238   if (strstr(verboseOption,"tim") || strstr(verboseOption,"all"))
 239     gBenchmark->Start("TOFSDigitizer");
 240
 241   if (fEdgeTails) ftail = new TF1("tail",TimeWithTail,-2,2,3);
 242
 243   Int_t nselectedHits=0;
 244   Int_t ntotalsdigits=0;
 245   Int_t ntotalupdates=0;
 246   Int_t nnoisesdigits=0;
 247   Int_t nsignalsdigits=0;
 248   Int_t nHitsFromPrim=0;
 249   Int_t nHitsFromSec=0;
 250   Int_t nlargeTofDiff=0;
 251
 252   Bool_t thereIsNotASelection=(fSelectedSector==-1) && (fSelectedPlate==-1);
 253
 254   if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
 255   gAlice = fRunLoader->GetAliRun();
 256
 257   fRunLoader->LoadKinematics();
 258
 259   AliTOF *tof = (AliTOF *) gAlice->GetDetector("TOF");
 260
 261   if (!tof) {
 262     AliError("TOF not found");
 263     return;
 264   }
 265
 266   fTOFLoader->LoadHits("read");
 267   fTOFLoader->LoadSDigits("recreate");
 268
 269   for (Int_t iEvent=fEvent1; iEvent<fEvent2; iEvent++) {
 270 //     cout << "------------------- "<< GetName() << " ------------- \n";
 271 //     cout << "Sdigitizing event " << iEvent << endl;
 272
 273     fRunLoader->GetEvent(iEvent);
 274
 275     TTree *hitTree = fTOFLoader->TreeH ();
 276     if (!hitTree) return;
 277
 278     if (fTOFLoader->TreeS () == 0) fTOFLoader->MakeTree ("S");
 279
 280     //Make branch for digits
 281     tof->MakeBranch("S");
 282
 283     // recreate TClonesArray fSDigits - for backward compatibility
 284     if (tof->SDigits() == 0) {
 285       tof->CreateSDigitsArray();
 286     } else {
 287       tof->RecreateSDigitsArray();
 288     }
 289
 290     tof->SetTreeAddress();
 291
 292     Int_t version=tof->IsVersion();
 293
 294     Int_t nselectedHitsinEv=0;
 295     Int_t ntotalsdigitsinEv=0;
 296     Int_t ntotalupdatesinEv=0;
 297     Int_t nnoisesdigitsinEv=0;
 298     Int_t nsignalsdigitsinEv=0;
 299
 300     TParticle *particle;
 301     //AliTOFhit *tofHit;
 302     TClonesArray *tofHitArray = tof->Hits();
 303
 304     // create hit map
 305     AliTOFHitMap *hitMap = new AliTOFHitMap(tof->SDigits(), fTOFGeometry);
 306
 307     TBranch * tofHitsBranch = hitTree->GetBranch("TOF");
 308
 309     Int_t ntracks = static_cast<Int_t>(hitTree->GetEntries());
 310     for (Int_t track = 0; track < ntracks; track++)
 311     {
 312       gAlice->ResetHits();
 313       tofHitsBranch->GetEvent(track);
 314       particle = gAlice->GetMCApp()->Particle(track);
 315       Int_t nhits = tofHitArray->GetEntriesFast();
 316       // cleaning all hits of the same track in the same pad volume
 317       // it is a rare event, however it happens
 318
 319       Int_t previousTrack =-1;
 320       Int_t previousSector=-1;
 321       Int_t previousPlate =-1;
 322       Int_t previousStrip =-1;
 323       Int_t previousPadX  =-1;
 324       Int_t previousPadZ  =-1;
 325
 326       for (Int_t hit = 0; hit < nhits; hit++) {
 327         Int_t    vol[5];       // location for a digit
 328         Float_t  digit[2];     // TOF digit variables
 329         Int_t tracknum;
 330         Float_t dxPad;
 331         Float_t dzPad;
 332         Float_t geantTime;
 333
 334         // fp: really sorry for this, it is a temporary trick to have
 335         // track length too
 336         if(version!=6 && version!=7){
 337           AliTOFhit *tofHit = (AliTOFhit *) tofHitArray->UncheckedAt(hit);
 338           tracknum = tofHit->GetTrack();
 339           vol[0] = tofHit->GetSector();
 340           vol[1] = tofHit->GetPlate();
 341           vol[2] = tofHit->GetStrip();
 342           vol[3] = tofHit->GetPadx();
 343           vol[4] = tofHit->GetPadz();
 344           dxPad = tofHit->GetDx();
 345           dzPad = tofHit->GetDz();
 346           geantTime = tofHit->GetTof(); // unit [s]
 347         } else {
 348           AliTOFhitT0 *tofHit = (AliTOFhitT0 *) tofHitArray->UncheckedAt(hit);
 349           tracknum = tofHit->GetTrack();
 350           vol[0] = tofHit->GetSector();
 351           vol[1] = tofHit->GetPlate();
 352           vol[2] = tofHit->GetStrip();
 353           vol[3] = tofHit->GetPadx();
 354           vol[4] = tofHit->GetPadz();
 355           dxPad = tofHit->GetDx();
 356           dzPad = tofHit->GetDz();
 357           geantTime = tofHit->GetTof(); // unit [s]
 358         }
 359
 360         geantTime *= 1.e+09;  // conversion from [s] to [ns]
 361
 362         // selection case for sdigitizing only hits in a given plate of a given sector
 363         if(thereIsNotASelection || (vol[0]==fSelectedSector && vol[1]==fSelectedPlate)){
 364
 365           Bool_t dummy=((tracknum==previousTrack) && (vol[0]==previousSector) && (vol[1]==previousPlate) && (vol[2]==previousStrip));
 366
 367           Bool_t isCloneOfThePrevious=dummy && ((vol[3]==previousPadX) && (vol[4]==previousPadZ));
 368
 369           Bool_t isNeighOfThePrevious=dummy && ((((vol[3]==previousPadX-1) || (vol[3]==previousPadX+1)) && (vol[4]==previousPadZ)) || ((vol[3]==previousPadX) && ((vol[4]==previousPadZ+1) || (vol[4]==previousPadZ-1))));
 370
 371           if(!isCloneOfThePrevious && !isNeighOfThePrevious){
 372             // update "previous" values
 373             // in fact, we are yet in the future, so the present is past
 374             previousTrack=tracknum;
 375             previousSector=vol[0];
 376             previousPlate=vol[1];
 377             previousStrip=vol[2];
 378             previousPadX=vol[3];
 379             previousPadZ=vol[4];
 380
 381             nselectedHits++;
 382             nselectedHitsinEv++;
 383             if (particle->GetFirstMother() < 0) nHitsFromPrim++; // counts hits due to primary particles
 384
 385             Float_t xStrip=AliTOFGeometry::XPad()*(vol[3]+0.5-0.5*AliTOFGeometry::NpadX())+dxPad;
 386             Float_t zStrip=AliTOFGeometry::ZPad()*(vol[4]+0.5-0.5*AliTOFGeometry::NpadZ())+dzPad;
 387
 388             Int_t nActivatedPads = 0, nFiredPads = 0;
 389             Bool_t isFired[4] = {kFALSE, kFALSE, kFALSE, kFALSE};
 390             Float_t tofAfterSimul[4] = {0., 0., 0., 0.};
 391             Float_t qInduced[4] = {0.,0.,0.,0.};
 392             Int_t nPlace[4] = {0, 0, 0, 0};
 393             Float_t averageTime = 0.;
 394             SimulateDetectorResponse(zStrip,xStrip,geantTime,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime);
 395             if(nFiredPads) {
 396               for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
 397                 if(isFired[indexOfPad]){ // the pad has fired
 398                   Float_t timediff=geantTime-tofAfterSimul[indexOfPad];
 399
 400                   if(timediff>=0.2) nlargeTofDiff++;
 401
 402                   digit[0] = (Int_t) ((tofAfterSimul[indexOfPad]*1.e+03)/AliTOFGeometry::TdcBinWidth()); // TDC bin number (each bin -> 24.4 ps)
 403
 404                   Float_t landauFactor = gRandom->Landau(fAdcMean, fAdcRms);
 405                   digit[1] = (Int_t) (qInduced[indexOfPad] * landauFactor); // ADC bins (each bin -> 0.25 (or 0.03) pC)
 406
 407                   // recalculate the volume only for neighbouring pads
 408                   if(indexOfPad){
 409                     (nPlace[indexOfPad]<=AliTOFGeometry::NpadX()) ? vol[4] = 0 : vol[4] = 1;
 410                     (nPlace[indexOfPad]<=AliTOFGeometry::NpadX()) ? vol[3] = nPlace[indexOfPad] - 1 : vol[3] = nPlace[indexOfPad] - AliTOFGeometry::NpadX() - 1;
 411                   }
 412                   // check if two sdigit are on the same pad;
 413                   // in that case we sum the two or more sdigits
 414                   if (hitMap->TestHit(vol) != kEmpty) {
 415                     AliTOFSDigit *sdig = static_cast<AliTOFSDigit*>(hitMap->GetHit(vol));
 416                     Int_t tdctime = (Int_t) digit[0];
 417                     Int_t adccharge = (Int_t) digit[1];
 418                     sdig->Update(AliTOFGeometry::TdcBinWidth(),tdctime,adccharge,tracknum);
 419                     ntotalupdatesinEv++;
 420                     ntotalupdates++;
 421                   } else {
 422
 423                     tof->AddSDigit(tracknum, vol, digit);
 424
 425                     if(indexOfPad){
 426                       nnoisesdigits++;
 427                       nnoisesdigitsinEv++;
 428                     } else {
 429                       nsignalsdigits++;
 430                       nsignalsdigitsinEv++;
 431                     }
 432                     ntotalsdigitsinEv++;
 433                     ntotalsdigits++;
 434                     hitMap->SetHit(vol);
 435                   } // if (hitMap->TestHit(vol) != kEmpty)
 436                 } // if(isFired[indexOfPad])
 437               } // end loop on nActivatedPads
 438             } // if(nFiredPads) i.e. if some pads has fired
 439           } // close if(!isCloneOfThePrevious)
 440         } // close the selection on sector and plate
 441       } // end loop on hits for the current track
 442     } // end loop on ntracks
 443
 444     delete hitMap;
 445
 446     fTOFLoader->TreeS()->Reset();
 447     fTOFLoader->TreeS()->Fill();
 448     fTOFLoader->WriteSDigits("OVERWRITE");
 449
 450     if (tof->SDigits()) tof->ResetSDigits();
 451
 452     if (strstr(verboseOption,"all")) {
 453       AliInfo("----------------------------------------");
 454       AliInfo("       <AliTOFSDigitizer>    ");
 455       AliInfo(Form("After sdigitizing %d hits in event %d", nselectedHitsinEv, iEvent));
 456       //" (" << nHitsFromPrim << " from primaries and " << nHitsFromSec << " from secondaries) TOF hits, "
 457       AliInfo(Form("%d digits have been created", ntotalsdigitsinEv));
 458       AliInfo(Form("(%d due to signals and %d due to border effect)", nsignalsdigitsinEv, nnoisesdigitsinEv));
 459       AliInfo(Form("%d total updates of the hit map have been performed in current event", ntotalupdatesinEv));
 460       AliInfo("----------------------------------------");
 461     }
 462
 463   } //event loop on events
 464
 465     fTOFLoader->UnloadSDigits();
 466     fTOFLoader->UnloadHits();
 467     fRunLoader->UnloadKinematics();
 468     //fRunLoader->UnloadgAlice();
 469
 470   // free used memory
 471   if (ftail){
 472     delete ftail;
 473     ftail = 0;
 474   }
 475
 476   nHitsFromSec=nselectedHits-nHitsFromPrim;
 477   if(strstr(verboseOption,"all")){
 478     AliInfo("----------------------------------------");
 479     AliInfo("----------------------------------------");
 480     AliInfo("-----------SDigitization Summary--------");
 481     AliInfo("       <AliTOFSDigitizer>     ");
 482     AliInfo(Form("After sdigitizing %d hits", nselectedHits));
 483     AliInfo(Form("in %d events", fEvent2-fEvent1));
 484 //" (" << nHitsFromPrim << " from primaries and " << nHitsFromSec << " from secondaries) TOF hits, "
 485     AliInfo(Form("%d sdigits have been created", ntotalsdigits));
 486     AliInfo(Form("(%d due to signals and "
 487                  "%d due to border effect)", nsignalsdigits, nnoisesdigits));
 488     AliInfo(Form("%d total updates of the hit map have been performed", ntotalupdates));
 489     AliInfo(Form("in %d cases the time of flight difference is greater than 200 ps", nlargeTofDiff));
 490   }
 491
 492
 493   if(strstr(verboseOption,"tim") || strstr(verboseOption,"all")){
 494     gBenchmark->Stop("TOFSDigitizer");
 495     AliInfo("AliTOFSDigitizer:");
 496     AliInfo(Form("   took %f seconds in order to make sdigits "
 497          "%f seconds per event", gBenchmark->GetCpuTime("TOFSDigitizer"), gBenchmark->GetCpuTime("TOFSDigitizer")/(fEvent2-fEvent1)));
 498     AliInfo(" +++++++++++++++++++++++++++++++++++++++++++++++++++ ");
 499   }
 500
 501 }
 502
 503 //__________________________________________________________________
 504 void AliTOFSDigitizer::Print(Option_t* /*opt*/)const
 505 {
 506   cout << "------------------- "<< GetName() << " ------------- \n";
 507 }
 508
 509 //__________________________________________________________________
 510 void AliTOFSDigitizer::SelectSectorAndPlate(Int_t sector, Int_t plate)
 511 {
 512   //Select sector and plate
 513   Bool_t isaWrongSelection=(sector < 0) || (sector >= AliTOFGeometry::NSectors()) || (plate < 0) || (plate >= AliTOFGeometry::NPlates());
 514   if(isaWrongSelection){
 515     AliError("You have selected an invalid value for sector or plate ");
 516     AliError(Form("The correct range for sector is [0,%d]", AliTOFGeometry::NSectors()-1));
 517     AliError(Form("The correct range for plate  is [0,%d]",  AliTOFGeometry::NPlates()-1));
 518     AliError("By default we continue sdigitizing all hits in all plates of all sectors");
 519   } else {
 520     fSelectedSector=sector;
 521     fSelectedPlate =plate;
 522     AliInfo(Form("SDigitizing only hits in plate %d of the sector %d", fSelectedPlate, fSelectedSector));
 523   }
 524 }
 525
 526 //__________________________________________________________________
 527 void AliTOFSDigitizer::SimulateDetectorResponse(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)
 528 {
 529   // Description:
 530   // Input:  z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
 531   //         geantTime - time generated by Geant, ns
 532   // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
 533   //         nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
 534   //         qInduced[iPad]- charge induced on pad, arb. units
 535   //                         this array is initialized at zero by the caller
 536   //         tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
 537   //                                   this array is initialized at zero by the caller
 538   //         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.
 539   //                       The weight is given by the qInduced[iPad]/qCenterPad
 540   //                                   this variable is initialized at zero by the caller
 541   //         nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
 542   //                                   this variable is initialized at zero by the caller
 543   //
 544   // Description of used variables:
 545   //         eff[iPad] - efficiency of the pad
 546   //         res[iPad] - resolution of the pad, ns
 547   //         timeWalk[iPad] - time walk of the pad, ns
 548   //         timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
 549   //         PadId[iPad] - Pad Identifier
 550   //                    E | F    -->   PadId[iPad] = 5 | 6
 551   //                    A | B    -->   PadId[iPad] = 1 | 2
 552   //                    C | D    -->   PadId[iPad] = 3 | 4
 553   //         nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
 554   //         qCenterPad - charge extimated for each pad, arb. units
 555   //         weightsSum - sum of weights extimated for each pad fired, arb. units
 556
 557   const Float_t kSigmaForTail[2] = {AliTOFGeometry::SigmaForTail1(),AliTOFGeometry::SigmaForTail2()}; //for tail
 558   Int_t iz = 0, ix = 0;
 559   Float_t dX = 0., dZ = 0., x = 0., z = 0.;
 560   Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
 561   Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
 562   Float_t logOfqInd = 0.;
 563   Float_t weightsSum = 0.;
 564   Int_t nTail[4]  = {0,0,0,0};
 565   Int_t padId[4]  = {0,0,0,0};
 566   Float_t eff[4]  = {0.,0.,0.,0.};
 567   Float_t res[4]  = {0.,0.,0.,0.};
 568   //  Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
 569   Float_t qCenterPad = 1.;
 570   Float_t timeWalk[4]  = {0.,0.,0.,0.};
 571   Float_t timeDelay[4] = {0.,0.,0.,0.};
 572
 573   nActivatedPads = 0;
 574   nFiredPads = 0;
 575
 576   (z0 <= 0) ? iz = 0 : iz = 1;
 577   dZ = z0 + (0.5 * AliTOFGeometry::NpadZ() - iz - 0.5) * AliTOFGeometry::ZPad(); // hit position in the pad frame, (0,0) - center of the pad
 578   z = 0.5 * AliTOFGeometry::ZPad() - TMath::Abs(dZ);                               // variable for eff., res. and timeWalk. functions
 579   iz++;                                                                              // z row: 1, ..., AliTOFGeometry::NpadZ = 2
 580   ix = (Int_t)((x0 + 0.5 * AliTOFGeometry::NpadX() * AliTOFGeometry::XPad()) / AliTOFGeometry::XPad());
 581   dX = x0 + (0.5 * AliTOFGeometry::NpadX() - ix - 0.5) * AliTOFGeometry::XPad(); // hit position in the pad frame, (0,0) - center of the pad
 582   x = 0.5 * AliTOFGeometry::XPad() - TMath::Abs(dX);                               // variable for eff., res. and timeWalk. functions;
 583   ix++;                                                                              // x row: 1, ..., AliTOFGeometry::NpadX = 48
 584
 585   ////// Pad A:
 586   nActivatedPads++;
 587   nPlace[nActivatedPads-1] = (iz - 1) * AliTOFGeometry::NpadX() + ix;
 588   qInduced[nActivatedPads-1] = qCenterPad;
 589   padId[nActivatedPads-1] = 1;
 590
 591   if (fEdgeEffect == 0) {
 592     eff[nActivatedPads-1] = fEffCenter;
 593     if (gRandom->Rndm() < eff[nActivatedPads-1]) {
 594       nFiredPads = 1;
 595       res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + fResCenter * fResCenter); // ns
 596       isFired[nActivatedPads-1] = kTRUE;
 597       tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
 598       averageTime = tofTime[nActivatedPads-1];
 599     }
 600   } else {
 601
 602     if(z < h) {
 603       if(z < h2) {
 604         effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
 605       } else {
 606         effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
 607       }
 608       resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
 609       timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
 610       nTail[nActivatedPads-1] = 1;
 611     } else {
 612       effZ = fEffCenter;
 613       resZ = fResCenter;
 614       timeWalkZ = fTimeWalkCenter;
 615     }
 616
 617     if(x < h) {
 618       if(x < h2) {
 619         effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
 620       } else {
 621         effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
 622       }
 623       resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
 624       timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
 625       nTail[nActivatedPads-1] = 1;
 626     } else {
 627       effX = fEffCenter;
 628       resX = fResCenter;
 629       timeWalkX = fTimeWalkCenter;
 630     }
 631
 632     (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
 633     (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
 634     (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 *  timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 635
 636
 637     ////// Pad B:
 638     if(z < k2) {
 639       effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
 640     } else {
 641       effZ = fEff3Boundary * (k - z) / (k - k2);
 642     }
 643     resZ = fResBoundary + fResSlope * z / k;
 644     timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
 645
 646     if(z < k && z > 0) {
 647       if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
 648         nActivatedPads++;
 649         nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX();
 650         eff[nActivatedPads-1] = effZ;
 651         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
 652         timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
 653         nTail[nActivatedPads-1] = 2;
 654         if (fTimeDelayFlag) {
 655           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
 656           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
 657           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
 658           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
 659           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 660         } else {
 661           timeDelay[nActivatedPads-1] = 0.;
 662         }
 663         padId[nActivatedPads-1] = 2;
 664       }
 665     }
 666
 667
 668     ////// Pad C, D, E, F:
 669     if(x < k2) {
 670       effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
 671     } else {
 672       effX = fEff3Boundary * (k - x) / (k - k2);
 673     }
 674     resX = fResBoundary + fResSlope*x/k;
 675     timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
 676
 677     if(x < k && x > 0) {
 678       //   C:
 679       if(ix > 1 && dX < 0) {
 680         nActivatedPads++;
 681         nPlace[nActivatedPads-1] = nPlace[0] - 1;
 682         eff[nActivatedPads-1] = effX;
 683         res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX); // ns
 684         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 685         nTail[nActivatedPads-1] = 2;
 686         if (fTimeDelayFlag) {
 687           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
 688           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
 689           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
 690           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
 691           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 692         } else {
 693           timeDelay[nActivatedPads-1] = 0.;
 694         }
 695         padId[nActivatedPads-1] = 3;
 696
 697         //     D:
 698         if(z < k && z > 0) {
 699           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
 700             nActivatedPads++;
 701             nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() - 1;
 702             eff[nActivatedPads-1] = effX * effZ;
 703             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
 704             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 705
 706             nTail[nActivatedPads-1] = 2;
 707             if (fTimeDelayFlag) {
 708               if (TMath::Abs(x) < TMath::Abs(z)) {
 709                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
 710                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
 711                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
 712                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
 713               } else {
 714                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
 715                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
 716                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
 717                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
 718               }
 719               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 720             } else {
 721               timeDelay[nActivatedPads-1] = 0.;
 722             }
 723             padId[nActivatedPads-1] = 4;
 724           }
 725         }  // end D
 726       }  // end C
 727
 728       //   E:
 729       if(ix < AliTOFGeometry::NpadX() && dX > 0) {
 730         nActivatedPads++;
 731         nPlace[nActivatedPads-1] = nPlace[0] + 1;
 732         eff[nActivatedPads-1] = effX;
 733         res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(fAddTRes*fAddTRes + resX * resX)); // ns
 734         timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
 735         nTail[nActivatedPads-1] = 2;
 736         if (fTimeDelayFlag) {
 737           //      qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
 738           //      qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
 739           qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
 740           logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
 741           timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 742         } else {
 743           timeDelay[nActivatedPads-1] = 0.;
 744         }
 745         padId[nActivatedPads-1] = 5;
 746
 747
 748         //     F:
 749         if(z < k && z > 0) {
 750           if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
 751             nActivatedPads++;
 752             nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFGeometry::NpadX() + 1;
 753             eff[nActivatedPads - 1] = effX * effZ;
 754             (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(fAddTRes*fAddTRes + resZ * resZ); // ns
 755             (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
 756             nTail[nActivatedPads-1] = 2;
 757             if (fTimeDelayFlag) {
 758               if (TMath::Abs(x) < TMath::Abs(z)) {
 759                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
 760                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
 761                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
 762                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
 763               } else {
 764                 //              qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
 765                 //              qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
 766                 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
 767                 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
 768               }
 769               timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
 770             } else {
 771               timeDelay[nActivatedPads-1] = 0.;
 772             }
 773             padId[nActivatedPads-1] = 6;
 774           }
 775         }  // end F
 776       }  // end E
 777     } // end if(x < k)
 778
 779
 780     for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
 781       if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
 782       if(gRandom->Rndm() < eff[iPad]) {
 783         isFired[iPad] = kTRUE;
 784         nFiredPads++;
 785         if(fEdgeTails) {
 786           if(nTail[iPad] == 0) {
 787             tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
 788           } else {
 789             ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
 790             Double_t timeAB = ftail->GetRandom();
 791             tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
 792           }
 793         } else {
 794           tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
 795         }
 796         if (fAverageTimeFlag) {
 797           averageTime += tofTime[iPad] * qInduced[iPad];
 798           weightsSum += qInduced[iPad];
 799         } else {
 800           averageTime += tofTime[iPad];
 801           weightsSum += 1.;
 802         }
 803       }
 804     }
 805     if (weightsSum!=0) averageTime /= weightsSum;
 806   } // end else (fEdgeEffect != 0)
 807 }
 808
 809 //__________________________________________________________________
 810 void AliTOFSDigitizer::PrintParameters()const
 811 {
 812   //
 813   // Print parameters used for sdigitization
 814   //
 815   cout << " ------------------- "<< GetName() << " -------------" << endl ;
 816   cout << " Parameters used for TOF SDigitization " << endl ;
 817   //  Printing the parameters
 818
 819   cout << " Number of events:                        " << (fEvent2-fEvent1) << endl;
 820   cout << " from event " << fEvent1 << " to event " << (fEvent2-1) << endl;
 821   cout << " Time Resolution (ns) "<< fTimeResolution <<" Pad Efficiency: "<< fpadefficiency << endl;
 822   cout << " Edge Effect option:  "<<  fEdgeEffect<< endl;
 823
 824   cout << " Boundary Effect Simulation Parameters " << endl;
 825   cout << " Hparameter: "<< fHparameter<<"  H2parameter:"<< fH2parameter <<"  Kparameter:"<< fKparameter<<"  K2parameter: "<< fK2parameter << endl;
 826   cout << " Efficiency in the central region of the pad: "<< fEffCenter << endl;
 827   cout << " Efficiency at the boundary region of the pad: "<< fEffBoundary << endl;
 828   cout << " Efficiency value at H2parameter "<< fEff2Boundary << endl;
 829   cout << " Efficiency value at K2parameter "<< fEff3Boundary << endl;
 830   cout << " Resolution (ps) in the central region of the pad: "<< fResCenter << endl;
 831   cout << " Resolution (ps) at the boundary of the pad      : "<< fResBoundary << endl;
 832   cout << " Slope (ps/K) for neighbouring pad               : "<< fResSlope <<endl;
 833   cout << " Time walk (ps) in the central region of the pad : "<< fTimeWalkCenter << endl;
 834   cout << " Time walk (ps) at the boundary of the pad       : "<< fTimeWalkBoundary<< endl;
 835   cout << " Slope (ps/K) for neighbouring pad               : "<< fTimeWalkSlope<<endl;
 836   cout << " Pulse Heigth Simulation Parameters " << endl;
 837   cout << " Flag for delay due to the PulseHeightEffect  : "<< fTimeDelayFlag <<endl;
 838   cout << " Pulse Height Slope                           : "<< fPulseHeightSlope<<endl;
 839   cout << " Time Delay Slope                             : "<< fTimeDelaySlope<<endl;
 840   cout << " Minimum charge amount which could be induced : "<< fMinimumCharge<<endl;
 841   cout << " Smearing in charge in (q1/q2) vs x plot      : "<< fChargeSmearing<<endl;
 842   cout << " Smearing in log of charge ratio              : "<< fLogChargeSmearing<<endl;
 843   cout << " Smearing in time in time vs log(q1/q2) plot  : "<< fTimeSmearing<<endl;
 844   cout << " Flag for average time                        : "<< fAverageTimeFlag<<endl;
 845   cout << " Edge tails option                            : "<< fEdgeTails << endl;
 846
 847 }