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