ITS/AliITSsimulationSDD.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 #include <Riostream.h>
  19 #include <stdlib.h>
  20 #include <stdio.h>
  21 #include <cstring>
  22
  23 #include <TCanvas.h>
  24 #include <TF1.h>
  25 #include <TH1.h>
  26 #include <TFile.h>
  27 #include <TRandom.h>
  28 #include <TROOT.h>
  29 #include "AliITS.h"
  30 #include "AliITSMapA2.h"
  31 #include "AliITSRawData.h"
  32 #include "AliITSdigitSPD.h"
  33 #include "AliITSetfSDD.h"
  34 #include "AliITSmodule.h"
  35 #include "AliITShit.h"
  36 #include "AliITSpList.h"
  37 #include "AliITSCalibrationSDD.h"
  38 #include "AliITSresponseSDD.h"
  39 #include "AliITSsimulationSDD.h"
  40 #include "AliLog.h"
  41 #include "AliRun.h"
  42
  43 ClassImp(AliITSsimulationSDD)
  44 ////////////////////////////////////////////////////////////////////////
  45 // Version: 0                                                         //
  46 // Written by Piergiorgio Cerello                                     //
  47 // November 23 1999                                                   //
  48 //                                                                    //
  49 // AliITSsimulationSDD is the simulation of SDDs.                     //
  50 ////////////////////////////////////////////////////////////////////////
  51
  52 //______________________________________________________________________
  53 AliITSsimulationSDD::AliITSsimulationSDD():
  54 AliITSsimulation(),
  55 fITS(0),
  56 fHitMap2(0),
  57 fHitSigMap2(0),
  58 fHitNoiMap2(0),
  59 fElectronics(0),
  60 fInZR(0),
  61 fInZI(0),
  62 fOutZR(0),
  63 fOutZI(0),
  64 fAnodeFire(0),
  65 fHis(0),
  66 fFlag(kFALSE),
  67 fCrosstalkFlag(kFALSE),
  68 fDoFFT(1),
  69 fNofMaps(0),
  70 fMaxNofSamples(0),
  71 fScaleSize(0){
  72     // Default constructor
  73     SetPerpendTracksFlag();
  74     SetCrosstalkFlag();
  75     SetDoFFT();
  76 }
  77 //______________________________________________________________________
  78 AliITSsimulationSDD::AliITSsimulationSDD(const AliITSsimulationSDD &source) :
  79     AliITSsimulation(source),
  80 fITS(source.fITS),
  81 fHitMap2(source.fHitMap2),
  82 fHitSigMap2(source.fHitSigMap2),
  83 fHitNoiMap2(source.fHitNoiMap2),
  84 fElectronics(source.fElectronics),
  85 fInZR(source.fInZR),
  86 fInZI(source.fInZI),
  87 fOutZR(source.fOutZR),
  88 fOutZI(source.fOutZI),
  89 fAnodeFire(source.fAnodeFire),
  90 fHis(source.fHis),
  91 fFlag(source.fFlag),
  92 fCrosstalkFlag(source.fCrosstalkFlag),
  93 fDoFFT(source.fDoFFT),
  94 fNofMaps(source.fNofMaps),
  95 fMaxNofSamples(source.fMaxNofSamples),
  96 fScaleSize(source.fScaleSize){
  97     // Copy constructor to satify Coding roules only.
  98
  99 }
 100 //______________________________________________________________________
 101 AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
 102     // Assignment operator to satify Coding roules only.
 103
 104     if(this==&src) return *this;
 105     Error("AliITSsimulationSDD","Not allowed to make a = with "
 106           "AliITSsimulationSDD Using default creater instead");
 107     return *this ;
 108 }
 109 /*
 110 //______________________________________________________________________
 111 AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
 112     // Assignment operator to satify Coding roules only.
 113
 114     if(this==&src) return *this;
 115     Error("AliITSsimulationSSD","Not allowed to make a = with "
 116           "AliITSsimulationSDD Using default creater instead");
 117     return *this ;
 118 }
 119 */
 120 //______________________________________________________________________
 121 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
 122 AliITSsimulation(dettyp),
 123 fITS(0),
 124 fHitMap2(0),
 125 fHitSigMap2(0),
 126 fHitNoiMap2(0),
 127 fElectronics(0),
 128 fInZR(0),
 129 fInZI(0),
 130 fOutZR(0),
 131 fOutZI(0),
 132 fAnodeFire(0),
 133 fHis(0),
 134 fFlag(kFALSE),
 135 fCrosstalkFlag(kFALSE),
 136 fDoFFT(1),
 137 fNofMaps(0),
 138 fMaxNofSamples(0),
 139 fScaleSize(0){
 140     // Default Constructor
 141   Init();
 142 }
 143 //______________________________________________________________________
 144 void AliITSsimulationSDD::Init(){
 145     // Standard Constructor
 146
 147     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 148     fScaleSize = ScaleFourier(seg);
 149     SetPerpendTracksFlag();
 150     SetCrosstalkFlag();
 151     SetDoFFT();
 152
 153     AliITSSimuParam* simpar = fDetType->GetSimuParam();
 154     fpList = new AliITSpList( seg->Npz(),
 155                               fScaleSize*seg->Npx() );
 156     fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
 157     fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
 158     fHitMap2 = fHitSigMap2;
 159
 160     fNofMaps = seg->Npz();
 161     fMaxNofSamples = seg->Npx();
 162     fAnodeFire = new Bool_t [fNofMaps];
 163
 164     Float_t sddWidth  = seg->Dz();
 165     Float_t anodePitch = seg->Dpz(0);
 166     Double_t timeStep  = (Double_t)seg->Dpx(0);
 167
 168     if(anodePitch*(fNofMaps/2) > sddWidth) {
 169       AliWarning(Form("Too many anodes %d or too big pitch %f ",
 170                 fNofMaps/2,anodePitch));
 171     } // end if
 172
 173
 174     fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
 175                                     simpar->GetSDDElectronics());
 176
 177
 178     fITS       = (AliITS*)gAlice->GetModule("ITS");
 179
 180     fInZR  = new Double_t [fScaleSize*fMaxNofSamples];
 181     fInZI  = new Double_t [fScaleSize*fMaxNofSamples];
 182     fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
 183     fOutZI = new Double_t [fScaleSize*fMaxNofSamples];
 184 }
 185 //______________________________________________________________________
 186 AliITSsimulationSDD::~AliITSsimulationSDD() {
 187     // destructor
 188
 189     //    delete fpList;
 190     delete fHitSigMap2;
 191     delete fHitNoiMap2;
 192     delete fElectronics;
 193
 194     fITS = 0;
 195
 196     if (fHis) {
 197         fHis->Delete();
 198         delete fHis;
 199     } // end if fHis
 200     if(fInZR)  delete [] fInZR;
 201     if(fInZI)  delete [] fInZI;
 202     if(fOutZR) delete [] fOutZR;
 203     if(fOutZI) delete [] fOutZI;
 204     if(fAnodeFire) delete [] fAnodeFire;
 205 }
 206 //______________________________________________________________________
 207 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
 208     // create maps to build the lists of tracks for each summable digit
 209     fModule = module;
 210     fEvent  = event;
 211     ClearMaps();
 212     memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);
 213 }
 214 //______________________________________________________________________
 215 void AliITSsimulationSDD::ClearMaps() {
 216     // clear maps
 217     fpList->ClearMap();
 218     fHitSigMap2->ClearMap();
 219     fHitNoiMap2->ClearMap();
 220 }
 221 //______________________________________________________________________
 222 void AliITSsimulationSDD::FastFourierTransform(Double_t *real,
 223                           Double_t *imag,Int_t direction) {
 224     // Do a Fast Fourier Transform
 225
 226     Int_t samples = fElectronics->GetSamples();
 227     Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
 228     Int_t m1 = samples;
 229     Int_t m  = samples/2;
 230     Int_t m2 = samples/m1;
 231     Int_t i,j,k;
 232     for(i=1; i<=l; i++) {
 233         for(j=0; j<samples; j += m1) {
 234             Int_t p = 0;
 235             for(k=j; k<= j+m-1; k++) {
 236                 Double_t wsr = fElectronics->GetWeightReal(p);
 237                 Double_t wsi = fElectronics->GetWeightImag(p);
 238                 if(direction == -1) wsi = -wsi;
 239                 Double_t xr = *(real+k+m);
 240                 Double_t xi = *(imag+k+m);
 241                 *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi);
 242                 *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr);
 243                 *(real+k) += xr;
 244                 *(imag+k) += xi;
 245                 p += m2;
 246             } // end for k
 247         } // end for j
 248         m1 = m;
 249         m /= 2;
 250         m2 += m2;
 251     } // end for i
 252     for(j=0; j<samples; j++) {
 253         Int_t j1 = j;
 254         Int_t p = 0;
 255         Int_t i1;
 256         for(i1=1; i1<=l; i1++) {
 257             Int_t j2 = j1;
 258             j1 /= 2;
 259             p = p + p + j2 - j1 - j1;
 260         } // end for i1
 261         if(p >= j) {
 262             Double_t xr = *(real+j);
 263             Double_t xi = *(imag+j);
 264             *(real+j) = *(real+p);
 265             *(imag+j) = *(imag+p);
 266             *(real+p) = xr;
 267             *(imag+p) = xi;
 268         } // end if p>=j
 269     } // end for j
 270     if(direction == -1) {
 271         for(i=0; i<samples; i++) {
 272             *(real+i) /= samples;
 273             *(imag+i) /= samples;
 274         } // end for i
 275     } // end if direction == -1
 276     return;
 277 }
 278
 279 //______________________________________________________________________
 280 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 281     // digitize module using the "slow" detector simulator creating
 282     // summable digits.
 283
 284     TObjArray *fHits = mod->GetHits();
 285     Int_t nhits      = fHits->GetEntriesFast();
 286     if( !nhits ) return;
 287
 288     InitSimulationModule( md, ev );
 289     HitsToAnalogDigits( mod );  // fills fHitMap2 which is = fHitSigmap2
 290     ChargeToSignal( fModule,kFALSE,kTRUE ); // - Process signal adding gain without adding noise
 291     fHitMap2 = fHitNoiMap2;   // - Swap to noise map
 292     ChargeToSignal( fModule,kTRUE,kFALSE );  // - Process only noise
 293     fHitMap2 = fHitSigMap2;   // - Return to signal map
 294     WriteSDigits();
 295     ClearMaps();
 296 }
 297 //______________________________________________________________________
 298 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
 299                                                Int_t mask ) {
 300     // Add Summable digits to module maps.
 301     AliITSSimuParam* simpar = fDetType->GetSimuParam();
 302     Int_t    nItems = pItemArray->GetEntries();
 303     Double_t maxadc = simpar->GetSDDMaxAdc();
 304     Bool_t sig = kFALSE;
 305
 306     // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
 307     for( Int_t i=0; i<nItems; i++ ) {
 308         AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
 309         if( pItem->GetModule() != fModule ) {
 310             Error( "AliITSsimulationSDD","Error reading, SDigits module "
 311                    "%d != current module %d: exit",
 312                    pItem->GetModule(), fModule );
 313             return sig;
 314         } // end if
 315
 316         if(pItem->GetSignal()>0.0 ) sig = kTRUE;
 317
 318         fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
 319         AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
 320         Double_t sigAE = pItem2->GetSignalAfterElect();
 321         if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
 322         Int_t ia;
 323         Int_t it;
 324         fpList->GetMapIndex( pItem->GetIndex(), ia, it );
 325         fHitMap2->SetHit( ia, it, sigAE );
 326         fAnodeFire[ia] = kTRUE;
 327     }
 328     return sig;
 329 }
 330 //______________________________________________________________________
 331 void AliITSsimulationSDD::FinishSDigitiseModule() {
 332     // digitize module using the "slow" detector simulator from
 333     // the sum of summable digits.
 334     FinishDigits() ;
 335     ClearMaps();
 336 }
 337 //______________________________________________________________________
 338 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 339     // create maps to build the lists of tracks for each digit
 340
 341     TObjArray *fHits = mod->GetHits();
 342     Int_t nhits      = fHits->GetEntriesFast();
 343
 344     InitSimulationModule( md, ev );
 345     if( !nhits ) return;
 346
 347     HitsToAnalogDigits( mod );
 348     ChargeToSignal( fModule,kTRUE,kTRUE );  // process signal + noise
 349
 350     for( Int_t i=0; i<fNofMaps; i++ ) {
 351         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 352             Int_t jdx = j*fScaleSize;
 353             Int_t index = fpList->GetHitIndex( i, j );
 354             AliITSpListItem pItemTmp2( fModule, index, 0. );
 355             // put the fScaleSize analog digits in only one
 356             for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 357                 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
 358                 if( pItemTmp == 0 ) continue;
 359                 pItemTmp2.Add( pItemTmp );
 360             }
 361             fpList->DeleteHit( i, j );
 362             fpList->AddItemTo( 0, &pItemTmp2 );
 363         }
 364     }
 365     FinishDigits();
 366     ClearMaps();
 367 }
 368 //______________________________________________________________________
 369 void AliITSsimulationSDD::FinishDigits() {
 370     // introduce the electronics effects and do zero-suppression if required
 371
 372     if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
 373
 374     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 375     Bool_t isZeroSupp = res->GetZeroSupp();
 376     if (isZeroSupp) Compress2D();
 377     else StoreAllDigits();
 378 }
 379 //______________________________________________________________________
 380 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
 381     // create maps to build the lists of tracks for each digit
 382   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 383   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 384   AliITSSimuParam* simpar = fDetType->GetSimuParam();
 385   TObjArray *hits     = mod->GetHits();
 386   Int_t      nhits    = hits->GetEntriesFast();
 387
 388   //    Int_t      arg[6]   = {0,0,0,0,0,0};
 389   Int_t     nofAnodes  = fNofMaps/2;
 390   Double_t  sddLength  = seg->Dx();
 391   Double_t  anodePitch = seg->Dpz(0);
 392   Double_t  timeStep   = seg->Dpx(0);
 393   Double_t  driftSpeed ;  // drift velocity (anode dependent)
 394   Double_t  nanoampToADC       = simpar->GetSDDMaxAdc()/simpar->GetSDDDynamicRange(); //   maxadc/topValue;
 395   Double_t  cHloss     = simpar->GetSDDChargeLoss();
 396   Float_t   dfCoeff, s1;
 397   simpar->GetSDDDiffCoeff(dfCoeff,s1); // Signal 2d Shape
 398   Double_t  eVpairs    = simpar->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
 399   Double_t  nsigma     = simpar->GetNSigmaIntegration(); //
 400   Int_t     nlookups   = simpar->GetGausNLookUp();       //
 401   Float_t   jitter     = simpar->GetSDDJitterError(); //
 402   Float_t   trigDelay  = simpar->GetSDDTrigDelay(); // compensation for MC time zero
 403   Float_t   timeZero=fDetType->GetResponseSDD()->GetTimeZero(fModule);
 404
 405   // Piergiorgio's part (apart for few variables which I made float
 406   // when i thought that can be done
 407   // Fill detector maps with GEANT hits
 408   // loop over hits in the module
 409
 410   const Float_t kconv = 1.0e+6;  // GeV->KeV
 411   Int_t     itrack      = 0;
 412   Int_t     iWing;       // which detector wing/side.
 413   Int_t     ii,kk,ka,kt; // loop indexs
 414   Int_t     ia,it,index; // sub-pixel integration indexies
 415   Int_t     iAnode;      // anode number.
 416   Int_t     timeSample;  // time buckett.
 417   Int_t     anodeWindow; // anode direction charge integration width
 418   Int_t     timeWindow;  // time direction charge integration width
 419   Int_t     jamin,jamax; // anode charge integration window
 420   Int_t     jtmin,jtmax; // time charge integration window
 421   Int_t     nsplitAn;    // the number of splits in anode and time windows
 422   Int_t     nsplitTb;    // the number of splits in anode and time windows
 423   Int_t     nOfSplits;   // number of times track length is split into
 424   Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
 425   Float_t   kkF;         // Floating point version of loop index kk.
 426   Double_t  pathInSDD; // Track length in SDD.
 427   Double_t  drPath; // average position of track in detector. in microns
 428   Double_t  drTime; // Drift time
 429   Double_t  avDrft;  // x position of path length segment in cm.
 430   Double_t  avAnode; // Anode for path length segment in Anode number (float)
 431   Double_t  zAnode;  // Floating point anode number.
 432   Double_t  driftPath; // avDrft in microns.
 433   Double_t  width;     // width of signal at anodes.
 434   Double_t  depEnergy; // Energy deposited in this GEANT step.
 435   Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
 436   Double_t  sigA; // sigma of signal at anode.
 437   Double_t  sigT; // sigma in time/drift direction for track segment
 438   Double_t  aStep,aConst; // sub-pixel size and offset anode
 439   Double_t  tStep,tConst; // sub-pixel size and offset time
 440   Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
 441   Double_t  chargeloss; // charge loss for track segment.
 442   Double_t  anodeAmplitude; // signal amplitude in anode direction
 443   Double_t  aExpo;          // exponent of Gaussian anode direction
 444   Double_t  timeAmplitude;  // signal amplitude in time direction
 445   Double_t  tExpo;          // exponent of Gaussian time direction
 446   Double_t  tof;            // Time of flight in ns of this step.
 447
 448   for(ii=0; ii<nhits; ii++) {
 449     if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
 450                           depEnergy,itrack)) continue;
 451     Float_t xloc=xL[0];
 452     if(xloc>0) iWing=0; // left side, carlos channel 0
 453     else iWing=1; // right side
 454
 455     Float_t zloc=xL[2]+0.5*dxL[2];
 456     zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
 457     driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
 458     if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
 459       AliWarning("Time Interval > Allowed Time Interval");
 460     }
 461     depEnergy  *= kconv;
 462     if (!depEnergy) {
 463       AliDebug(1,
 464                Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
 465                     itrack,ii,mod->GetIndex()));
 466       continue;
 467       // continue if the particle did not lose energy
 468       // passing through detector
 469     } // end if !depEnergy
 470
 471     tof=0.;
 472     AliITShit* h=(AliITShit*)hits->At(ii);
 473     if(h) tof=h->GetTOF()*1E9;
 474     AliDebug(1,Form("TOF for hit %d on mod %d (particle %d)=%g",ii,fModule,h->Track(),tof));
 475
 476     xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
 477     pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
 478
 479     if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
 480     drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
 481     drPath = sddLength-drPath;
 482     if(drPath < 0) {
 483       AliDebug(1, // this should be fixed at geometry level
 484                Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
 485                     drPath,sddLength,dxL[0],xL[0]));
 486       continue;
 487     } // end if drPath < 0
 488
 489     // Compute number of segments to brake step path into
 490     drTime = drPath/driftSpeed;  //   Drift Time
 491     sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
 492     // calcuate the number of time the path length should be split into.
 493     nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
 494     if(fFlag) nOfSplits = 1;
 495
 496     // loop over path segments, init. some variables.
 497     depEnergy /= nOfSplits;
 498     nOfSplitsF = (Float_t) nOfSplits;
 499     Float_t theAverage=0.,theSteps=0.;
 500     for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
 501       kkF       = (Float_t) kk + 0.5;
 502       avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
 503       avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
 504       theSteps+=1.;
 505       theAverage+=avAnode;
 506       zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
 507       driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
 508       driftPath = TMath::Abs(10000.*avDrft);
 509       driftPath = sddLength-driftPath;
 510       if(driftPath < 0) {
 511         AliDebug(1, // this should be fixed at geometry level
 512                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
 513                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
 514         continue;
 515       } // end if driftPath < 0
 516       drTime     = driftPath/driftSpeed; // drift time for segment.
 517       // Sigma along the anodes for track segment.
 518       sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
 519       sigT       = sigA/driftSpeed;
 520
 521       drTime+=tof; // take into account Time Of Flight from production point
 522       drTime-=trigDelay;
 523       drTime+=timeZero;
 524       timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1.001); // time bin in range 1-256 !!!
 525       if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
 526       iAnode = (Int_t) (1.001+zAnode); // iAnode in range 1-256 !!!!
 527
 528         // Peak amplitude in nanoAmpere
 529       amplitude  = fScaleSize*160.*depEnergy/
 530         (timeStep*eVpairs*2.*acos(-1.));
 531       chargeloss = 1.-cHloss*driftPath/1000.;
 532       amplitude *= chargeloss;
 533       width  = 2.*nsigma/(nlookups-1);
 534       // Spread the charge
 535       nsplitAn = 4;
 536       nsplitTb=4;
 537       aStep  = anodePitch/(nsplitAn*sigA);
 538       aConst = zAnode*anodePitch/sigA;
 539       tStep  = timeStep/(nsplitTb*fScaleSize*sigT);
 540       tConst = drTime/sigT;
 541       // Define SDD window corresponding to the hit
 542       anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
 543       timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
 544       jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
 545       if(jamin <= 0) jamin = 1;
 546       if(jamin > nofAnodes*nsplitAn){
 547         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode min=%d",jamin));
 548         continue;
 549       }
 550       jamax = (iAnode + anodeWindow + 2)*nsplitAn;
 551       if(jamax > nofAnodes*nsplitAn) jamax = nofAnodes*nsplitAn;
 552       if(jamax <=0){
 553         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode max=%d",jamax));
 554         continue;
 555       }
 556       jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1;
 557       if(jtmin <= 0) jtmin = 1;
 558       if(jtmin > fScaleSize*fMaxNofSamples*nsplitTb){
 559         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample min=%d  tof=%f",jtmin,tof));
 560         continue;
 561       }
 562       jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb;
 563       if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb) jtmax = fScaleSize*fMaxNofSamples*nsplitTb;
 564       if(jtmax <= 0){
 565         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample max=%d  tof=%f",jtmax,tof));
 566         continue;
 567       }
 568
 569       // Spread the charge in the anode-time window
 570       for(ka=jamin; ka <=jamax; ka++) {
 571         ia = (ka-1)/nsplitAn + 1;
 572         if(ia <= 0) ia=1;
 573         if(ia > nofAnodes) ia = nofAnodes;
 574         aExpo     = (aStep*(ka-0.5)-aConst);
 575         if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
 576         else {
 577           Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
 578           anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
 579         }
 580         // index starts from 0
 581         index = iWing*nofAnodes+ia-1;
 582         if(anodeAmplitude){
 583           for(kt=jtmin; kt<=jtmax; kt++) {
 584             it = (kt-1)/nsplitTb+1;  // it starts from 1
 585             if(it<=0) it=1;
 586             if(it>fScaleSize*fMaxNofSamples)
 587               it = fScaleSize*fMaxNofSamples;
 588             tExpo    = (tStep*(kt-0.5)-tConst);
 589             if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
 590             else {
 591               Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
 592               timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep;
 593             }
 594             timeAmplitude *= nanoampToADC;
 595             //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
 596             Double_t charge = timeAmplitude;
 597             charge += fHitMap2->GetSignal(index,it-1);
 598             fHitMap2->SetHit(index, it-1, charge);
 599             fpList->AddSignal(index,it-1,itrack,ii-1,
 600                               mod->GetIndex(),timeAmplitude);
 601             fAnodeFire[index] = kTRUE;
 602           }  // end loop over time in window
 603         } // end if anodeAmplitude
 604       } // loop over anodes in window
 605     } // end loop over "sub-hits"
 606   } // end loop over hits
 607 }
 608
 609 //____________________________________________
 610 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) {
 611   // Adds a Digit.
 612   Int_t size = AliITSdigit::GetNTracks();
 613
 614   Int_t digits[3];
 615   Int_t * tracks = new Int_t[size];
 616   Int_t * hits = new Int_t[size];
 617   Float_t phys;
 618   Float_t * charges = new Float_t[size];
 619
 620   digits[0] = i;
 621   digits[1] = j;
 622   digits[2] = signalc;
 623
 624   AliITSpListItem *pItem = fpList->GetpListItem( i, j );
 625   if( pItem == 0 ) {
 626     phys = 0.0;
 627     for( Int_t l=0; l<size; l++ ) {
 628       tracks[l]  = 0;
 629       hits[l]    = 0;
 630       charges[l] = 0.0;
 631     }
 632   } else {
 633     Int_t idtrack =  pItem->GetTrack( 0 );
 634     if( idtrack >= 0 ) phys = pItem->GetSignal();
 635     else phys = 0.0;
 636
 637     for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
 638       tracks[l]  = pItem->GetTrack( l );
 639       hits[l]    = pItem->GetHit( l );
 640       charges[l] = pItem->GetSignal( l );
 641     }else{
 642       tracks[l]  = -3;
 643       hits[l]    = -1;
 644       charges[l] = 0.0;
 645     }// end for if
 646   }
 647
 648   fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale );
 649   delete [] tracks;
 650   delete [] hits;
 651   delete [] charges;
 652 }
 653 //______________________________________________________________________
 654 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) {
 655   // add baseline, noise, gain, electronics and ADC saturation effects
 656   // apply dead channels
 657
 658   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 659   Double_t baseline=0;
 660   Double_t noise=0;
 661   Double_t gain=0;
 662   Float_t contrib=0;
 663   Int_t i,k,kk;
 664   AliITSSimuParam* simpar = fDetType->GetSimuParam();
 665   Float_t maxadc = simpar->GetSDDMaxAdc();
 666   Int_t nGroup=fScaleSize;
 667   if(res->IsAMAt20MHz()){
 668     nGroup=fScaleSize/2;
 669   }
 670
 671   for (i=0;i<fNofMaps;i++) {
 672     if( !fAnodeFire[i] ) continue;
 673     baseline = res->GetBaseline(i);
 674     noise = res->GetNoise(i);
 675     gain = res->GetChannelGain(i)/fDetType->GetAverageGainSDD();
 676     if(res->IsBad()) gain=0.;
 677     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
 678     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 679       fInZR[k]  = fHitMap2->GetSignal(i,k);
 680       if(bAddGain) fInZR[k]*=gain;
 681       if( bAddNoise ) {
 682         contrib   = (baseline + noise*gRandom->Gaus());
 683         fInZR[k] += contrib;
 684       }
 685       fInZI[k]  = 0.;
 686     } // end for k
 687     if(!fDoFFT) {
 688       for(k=0; k<fMaxNofSamples; k++) {
 689         Double_t newcont = 0.;
 690         Double_t maxcont = 0.;
 691         for(kk=0;kk<fScaleSize;kk++) {
 692           newcont = fInZR[fScaleSize*k+kk];
 693           if(newcont > maxcont) maxcont = newcont;
 694         } // end for kk
 695         newcont = maxcont;
 696         if (newcont >= maxadc) newcont = maxadc -1;
 697         if(newcont >= baseline){
 698           Warning("","newcont=%f>=baseline=%f",newcont,baseline);
 699         } // end if
 700           // back to analog: ?
 701         fHitMap2->SetHit(i,k,newcont);
 702       }  // end for k
 703     }else{
 704       FastFourierTransform(&fInZR[0],&fInZI[0],1);
 705       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 706         Double_t rw = fElectronics->GetTraFunReal(k);
 707         Double_t iw = fElectronics->GetTraFunImag(k);
 708         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 709         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 710       } // end for k
 711       FastFourierTransform(&fOutZR[0],&fOutZI[0],-1);
 712       for(k=0; k<fMaxNofSamples; k++) {
 713         Double_t newcont1 = 0.;
 714         Double_t maxcont1 = 0.;
 715         for(kk=0;kk<nGroup;kk++) {
 716           newcont1 = fOutZR[fScaleSize*k+kk];
 717           if(newcont1 > maxcont1) maxcont1 = newcont1;
 718         } // end for kk
 719         newcont1 = maxcont1;
 720         if (newcont1 >= maxadc) newcont1 = maxadc -1;
 721         fHitMap2->SetHit(i,k,newcont1);
 722       } // end for k
 723     }
 724   } // end for i loop over anodes
 725   return;
 726 }
 727
 728 //______________________________________________________________________
 729 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 730     // function add the crosstalk effect to signal
 731     // temporal function, should be checked...!!!
 732
 733     // create and inizialice crosstalk map
 734     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 735     if( ctk == NULL ) {
 736         Error( "ApplyCrosstalk", "no memory for temporal map: exit " );
 737         return;
 738     }
 739     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 740     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 741     for( Int_t z=0; z<fNofMaps; z++ ) {
 742       Double_t baseline = calibr->GetBaseline(z);
 743         Bool_t on = kFALSE;
 744         Int_t tstart = 0;
 745         Int_t tstop = 0;
 746         Int_t nTsteps = 0;
 747
 748         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 749             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 750             if( fadc > baseline ) {
 751                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 752                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 753                     if( fadc1 < fadc ) continue;
 754                     on = kTRUE;
 755                     nTsteps = 0;
 756                     tstart = l;
 757                 }
 758                 nTsteps++;
 759             }
 760             else { // end fadc > baseline
 761                 if( on == kTRUE ) {
 762                     if( nTsteps > 2 ) {
 763                         tstop = l;
 764                         // make smooth derivative
 765                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 766                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 767                         if( ctk == NULL ) {
 768                             Error( "ApplyCrosstalk",
 769                                    "no memory for temporal array: exit " );
 770                             return;
 771                         }
 772                         for( Int_t i=tstart; i<tstop; i++ ) {
 773                             if( i > 2 && i < fMaxNofSamples-2 )
 774                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
 775                                     -0.1*fHitMap2->GetSignal( z,i-1 )
 776                                     +0.1*fHitMap2->GetSignal( z,i+1 )
 777                                     +0.2*fHitMap2->GetSignal( z,i+2 );
 778                         }
 779
 780                         // add crosstalk contribution to neibourg anodes
 781                         for( Int_t i=tstart; i<tstop; i++ ) {
 782                             Int_t anode = z - 1;
 783                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
 784                             Float_t ctktmp =  -dev[i1] * 0.25;
 785                             if( anode > 0 ) {
 786                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 787                             }
 788                             anode = z + 1;
 789                             if( anode < fNofMaps ) {
 790                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 791                             }
 792                         }
 793                         delete [] dev;
 794
 795                     } // if( nTsteps > 2 )
 796                     on = kFALSE;
 797                 }  // if( on == kTRUE )
 798             }  // else
 799         }
 800     }
 801
 802     for( Int_t a=0; a<fNofMaps; a++ )
 803         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
 804             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
 805             fHitMap2->SetHit( a, t, signal );
 806         }
 807
 808     delete [] ctk;
 809 }
 810
 811 //______________________________________________________________________
 812 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
 813     // To the 10 to 8 bit lossive compression.
 814     // code from Davide C. and Albert W.
 815
 816     if (signal < 128)  return signal;
 817     if (signal < 256)  return (128+((signal-128)>>1));
 818     if (signal < 512)  return (192+((signal-256)>>3));
 819     if (signal < 1024) return (224+((signal-512)>>4));
 820     return 0;
 821 }
 822 //______________________________________________________________________
 823 Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const {
 824   // Decompression from 8 to 10 bit
 825
 826   if (signal < 0 || signal > 255) {
 827     AliWarning(Form("Signal value %d out of range",signal));
 828     return 0;
 829   } // end if signal <0 || signal >255
 830
 831   if (signal < 128) return signal;
 832   if (signal < 192) {
 833     if (TMath::Odd(signal)) return (128+((signal-128)<<1));
 834     else  return (128+((signal-128)<<1)+1);
 835   } // end if signal < 192
 836   if (signal < 224) {
 837     if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3);
 838     else  return (256+((signal-192)<<3)+4);
 839   } // end if signal < 224
 840   if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7);
 841   return (512+((signal-224)<<4)+8);
 842 }
 843 //______________________________________________________________________
 844 void AliITSsimulationSDD::Compress2D(){
 845   // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10
 846   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 847   for (Int_t iWing=0; iWing<2; iWing++) {
 848     Int_t tL=res->GetZSLowThreshold(iWing);
 849     Int_t tH=res->GetZSHighThreshold(iWing);
 850     for (Int_t i=0; i<fNofMaps/2; i++) {
 851       Int_t ian=i+iWing*fNofMaps/2;
 852       if( !fAnodeFire[ian] ) continue;
 853       for (Int_t itb=0; itb<fMaxNofSamples; itb++) {
 854         Int_t nLow=0, nHigh=0;
 855         Float_t cC=fHitMap2->GetSignal(ian,itb);
 856         if(cC<=tL) continue;
 857         nLow++; // cC is greater than tL
 858         if(cC>tH) nHigh++;
 859         //                     N
 860         // Get "quintuple":   WCE
 861         //                     S
 862         Float_t wW=0.;
 863         if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1);
 864         if(wW>tL) nLow++;
 865         if(wW>tH) nHigh++;
 866         Float_t eE=0.;
 867         if(itb<fMaxNofSamples-1) eE=fHitMap2->GetSignal(ian,itb+1);
 868         if(eE>tL) nLow++;
 869         if(eE>tH) nHigh++;
 870         Float_t nN=0.;
 871         if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb);
 872         if(nN>tL) nLow++;
 873         if(nN>tH) nHigh++;
 874         Float_t sS=0.;
 875         if(i>0) sS=fHitMap2->GetSignal(ian-1,itb);
 876         if(sS>tL) nLow++;
 877         if(sS>tH) nHigh++;
 878
 879         if(nLow>=2 && nHigh>=1){
 880           Int_t signal=(Int_t)cC;
 881           Int_t signalc = Convert10to8(signal);
 882           Int_t signale = Convert8to10(signalc);
 883           signalc-=tL; // subtract low threshold after 10 to 8 bit compression
 884           if(signalc>=4) AddDigit(ian,itb,signalc,signale);  // store C
 885         }
 886       }
 887     }
 888   }
 889 }
 890
 891
 892 //______________________________________________________________________
 893 void AliITSsimulationSDD::StoreAllDigits(){
 894   // store digits for non-zero-suppressed data
 895   for (Int_t ian=0; ian<fNofMaps; ian++) {
 896     for (Int_t itb=0; itb<fMaxNofSamples; itb++){
 897       Int_t signal=(Int_t)(fHitMap2->GetSignal(ian,itb));
 898       Int_t signalc = Convert10to8(signal);
 899       Int_t signale = Convert8to10(signalc);
 900       AddDigit(ian,itb,signalc,signale);
 901     }
 902   }
 903 }
 904 //______________________________________________________________________
 905 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
 906     // Creates histograms of maps for debugging
 907     Int_t i;
 908
 909     fHis=new TObjArray(fNofMaps);
 910     for (i=0;i<fNofMaps;i++) {
 911         TString sddName("sdd_");
 912         Char_t candNum[4];
 913         sprintf(candNum,"%d",i+1);
 914         sddName.Append(candNum);
 915         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
 916                              0.,(Float_t) scale*fMaxNofSamples), i);
 917     } // end for i
 918 }
 919 //______________________________________________________________________
 920 void AliITSsimulationSDD::FillHistograms(){
 921     // fill 1D histograms from map
 922
 923     if (!fHis) return;
 924
 925     for( Int_t i=0; i<fNofMaps; i++) {
 926         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
 927         Int_t nsamples = hist->GetNbinsX();
 928         for( Int_t j=0; j<nsamples; j++) {
 929             Double_t signal=fHitMap2->GetSignal(i,j);
 930             hist->Fill((Float_t)j,signal);
 931         } // end for j
 932     } // end for i
 933 }
 934 //______________________________________________________________________
 935 void AliITSsimulationSDD::ResetHistograms(){
 936     // Reset histograms for this detector
 937     Int_t i;
 938
 939     for (i=0;i<fNofMaps;i++ ) {
 940         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
 941     } // end for i
 942 }
 943 //______________________________________________________________________
 944 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
 945     // Fills a histogram from a give anode.
 946
 947     if (!fHis) return 0;
 948
 949     if(wing <=0 || wing > 2) {
 950         Warning("GetAnode","Wrong wing number: %d",wing);
 951         return NULL;
 952     } // end if wing <=0 || wing >2
 953     if(anode <=0 || anode > fNofMaps/2) {
 954         Warning("GetAnode","Wrong anode number: %d",anode);
 955         return NULL;
 956     } // end if ampde <=0 || andoe > fNofMaps/2
 957
 958     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
 959     return (TH1F*)(fHis->At(index));
 960 }
 961 //______________________________________________________________________
 962 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
 963     // Writes the histograms to a file
 964
 965     if (!fHis) return;
 966
 967     hfile->cd();
 968     Int_t i;
 969     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
 970     return;
 971 }
 972 //______________________________________________________________________
 973 void AliITSsimulationSDD::WriteSDigits(){
 974     // Fills the Summable digits Tree
 975     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 976
 977     for( Int_t i=0; i<fNofMaps; i++ ) {
 978         if( !fAnodeFire[i] ) continue;
 979         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 980             Double_t sig = fHitMap2->GetSignal( i, j );
 981             if( sig > 0.2 ) {
 982                 Int_t jdx = j*fScaleSize;
 983                 Int_t index = fpList->GetHitIndex( i, j );
 984                 AliITSpListItem pItemTmp2( fModule, index, 0. );
 985                 // put the fScaleSize analog digits in only one
 986                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 987                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
 988                     if( pItemTmp == 0 ) continue;
 989                     pItemTmp2.Add( pItemTmp );
 990                 }
 991                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
 992                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
 993                 aliITS->AddSumDigit( pItemTmp2 );
 994             } // end if (sig > 0.2)
 995         }
 996     }
 997     return;
 998 }
 999 //______________________________________________________________________
1000 void AliITSsimulationSDD::PrintStatus() const {
1001     // Print SDD simulation Parameters
1002
1003     cout << "**************************************************" << endl;
1004     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1005     cout << "**************************************************" << endl;
1006     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1007     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1008     cout << "Number of Anodes used: " << fNofMaps << endl;
1009     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1010     cout << "Scale size factor: " << fScaleSize << endl;
1011     cout << "**************************************************" << endl;
1012 }