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   mapsmear   = simpar->GetSDDCorrMapPrecision(); //
 403   Float_t   trigDelay  = simpar->GetSDDTrigDelay(); // compensation for MC time zero
 404   if(res->IsAMAt20MHz()) trigDelay+=12.5; // compensation for discretization step
 405
 406   Float_t   timeZero=fDetType->GetResponseSDD()->GetTimeZero(fModule);
 407   Float_t   adcscale   = fDetType->GetResponseSDD()->GetADCtokeV(fModule);
 408   adcscale/=simpar->GetSDDkeVtoADC();
 409
 410   // Piergiorgio's part (apart for few variables which I made float
 411   // when i thought that can be done
 412   // Fill detector maps with GEANT hits
 413   // loop over hits in the module
 414
 415   const Float_t kconv = 1.0e+6;  // GeV->KeV
 416   Int_t     itrack      = 0;
 417   Int_t     iWing;       // which detector wing/side.
 418   Int_t     ii,kk,ka,kt; // loop indexs
 419   Int_t     ia,it,index; // sub-pixel integration indexies
 420   Int_t     iAnode;      // anode number.
 421   Int_t     timeSample;  // time buckett.
 422   Int_t     anodeWindow; // anode direction charge integration width
 423   Int_t     timeWindow;  // time direction charge integration width
 424   Int_t     jamin,jamax; // anode charge integration window
 425   Int_t     jtmin,jtmax; // time charge integration window
 426   Int_t     nsplitAn;    // the number of splits in anode and time windows
 427   Int_t     nsplitTb;    // the number of splits in anode and time windows
 428   Int_t     nOfSplits;   // number of times track length is split into
 429   Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
 430   Float_t   kkF;         // Floating point version of loop index kk.
 431   Double_t  pathInSDD; // Track length in SDD.
 432   Double_t  drPath; // average position of track in detector. in microns
 433   Double_t  drTime; // Drift time
 434   Double_t  avDrft;  // x position of path length segment in cm.
 435   Double_t  avAnode; // Anode for path length segment in Anode number (float)
 436   Double_t  zAnode;  // Floating point anode number.
 437   Double_t  driftPath; // avDrft in microns.
 438   Double_t  width;     // width of signal at anodes.
 439   Double_t  depEnergy; // Energy deposited in this GEANT step.
 440   Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
 441   Double_t  sigA; // sigma of signal at anode.
 442   Double_t  sigT; // sigma in time/drift direction for track segment
 443   Double_t  aStep,aConst; // sub-pixel size and offset anode
 444   Double_t  tStep,tConst; // sub-pixel size and offset time
 445   Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
 446   Double_t  chargeloss; // charge loss for track segment.
 447   Double_t  anodeAmplitude; // signal amplitude in anode direction
 448   Double_t  aExpo;          // exponent of Gaussian anode direction
 449   Double_t  timeAmplitude;  // signal amplitude in time direction
 450   Double_t  tExpo;          // exponent of Gaussian time direction
 451   Double_t  tof;            // Time of flight in ns of this step.
 452
 453   for(ii=0; ii<nhits; ii++) {
 454     if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
 455                           depEnergy,itrack)) continue;
 456     Float_t xloc=xL[0];
 457     if(xloc>0) iWing=0; // left side, carlos channel 0
 458     else iWing=1; // right side
 459
 460     Float_t zloc=xL[2]+0.5*dxL[2];
 461     zAnode=seg->GetAnodeFromLocal(xloc,zloc); // anode number in the range 0.-511.
 462     driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
 463     driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
 464
 465     if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
 466       AliWarning("Time Interval > Allowed Time Interval");
 467     }
 468     depEnergy  *= kconv;
 469     if (!depEnergy) {
 470       AliDebug(1,
 471                Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
 472                     itrack,ii,mod->GetIndex()));
 473       continue;
 474       // continue if the particle did not lose energy
 475       // passing through detector
 476     } // end if !depEnergy
 477
 478     tof=0.;
 479     AliITShit* h=(AliITShit*)hits->At(ii);
 480     if(h){
 481       tof=h->GetTOF()*1E9;
 482       AliDebug(1,Form("TOF for hit %d on mod %d (particle %d)=%g",ii,fModule,h->Track(),tof));
 483     }
 484
 485     Float_t corrx=0, corrz=0;
 486     res->GetShiftsForSimulation(xL[2],xL[0],corrz,corrx,seg);
 487     xL[2]-=corrz;
 488     xL[0]-=corrx;
 489     xL[0] += 0.0001*gRandom->Gaus( 0, mapsmear); //
 490     xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
 491
 492     pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
 493
 494     if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
 495     drPath = TMath::Abs(10000.*(dxL[0]+2.*xL[0])*0.5);
 496     drPath = sddLength-drPath;
 497     if(drPath < 0) {
 498       AliInfo( // this should be fixed at geometry level
 499                Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
 500                     drPath,sddLength,dxL[0],xL[0]));
 501       continue;
 502     } // end if drPath < 0
 503
 504     // Compute number of segments to brake step path into
 505     drTime = drPath/driftSpeed;  //   Drift Time
 506     sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
 507     // calcuate the number of time the path length should be split into.
 508     nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
 509     if(fFlag) nOfSplits = 1;
 510
 511     // loop over path segments, init. some variables.
 512     depEnergy /= nOfSplits;
 513     nOfSplitsF = (Float_t) nOfSplits;
 514     Float_t theAverage=0.,theSteps=0.;
 515     for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
 516       kkF       = (Float_t) kk + 0.5;
 517       avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
 518       avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
 519       theSteps+=1.;
 520       theAverage+=avAnode;
 521       zAnode = seg->GetAnodeFromLocal(avDrft,avAnode);
 522       driftSpeed = res->GetDriftSpeedAtAnode(zAnode);
 523       driftSpeed+= fDetType->GetResponseSDD()->GetDeltaVDrift(fModule,zAnode>255);
 524       driftPath = TMath::Abs(10000.*avDrft);
 525       driftPath = sddLength-driftPath;
 526       if(driftPath < 0) {
 527         AliDebug(1, // this should be fixed at geometry level
 528                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
 529                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
 530         continue;
 531       } // end if driftPath < 0
 532       drTime     = driftPath/driftSpeed; // drift time for segment.
 533       // Sigma along the anodes for track segment.
 534       sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
 535       sigT       = sigA/driftSpeed;
 536
 537       drTime+=tof; // take into account Time Of Flight from production point
 538       drTime-=trigDelay;
 539       drTime+=timeZero;
 540       timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1.001); // time bin in range 1-256 !!!
 541       if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
 542       iAnode = (Int_t) (1.001+zAnode); // iAnode in range 1-256 !!!!
 543
 544         // Peak amplitude in nanoAmpere
 545       amplitude  = fScaleSize*160.*depEnergy/
 546         (timeStep*eVpairs*2.*acos(-1.));
 547       chargeloss = 1.-cHloss*driftPath/1000.;
 548       amplitude *= chargeloss;
 549       amplitude *= adcscale;
 550       width  = 2.*nsigma/(nlookups-1);
 551       // Spread the charge
 552       nsplitAn = 4;
 553       nsplitTb=4;
 554       aStep  = anodePitch/(nsplitAn*sigA);
 555       aConst = zAnode*anodePitch/sigA;
 556       tStep  = timeStep/(nsplitTb*fScaleSize*sigT);
 557       tConst = drTime/sigT;
 558       // Define SDD window corresponding to the hit
 559       anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
 560       timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
 561       jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
 562       if(jamin <= 0) jamin = 1;
 563       if(jamin > nofAnodes*nsplitAn){
 564         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode min=%d",jamin));
 565         continue;
 566       }
 567       jamax = (iAnode + anodeWindow + 2)*nsplitAn;
 568       if(jamax > nofAnodes*nsplitAn) jamax = nofAnodes*nsplitAn;
 569       if(jamax <=0){
 570         AliDebug(1,Form("Energy deposition completely outside anode acceptance: anode max=%d",jamax));
 571         continue;
 572       }
 573       jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1;
 574       if(jtmin <= 0) jtmin = 1;
 575       if(jtmin > fScaleSize*fMaxNofSamples*nsplitTb){
 576         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample min=%d  tof=%f",jtmin,tof));
 577         continue;
 578       }
 579       jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb;
 580       if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb) jtmax = fScaleSize*fMaxNofSamples*nsplitTb;
 581       if(jtmax <= 0){
 582         AliDebug(1,Form("Energy deposition completely outside time acceptance: time sample max=%d  tof=%f",jtmax,tof));
 583         continue;
 584       }
 585
 586       // Spread the charge in the anode-time window
 587       for(ka=jamin; ka <=jamax; ka++) {
 588         ia = (ka-1)/nsplitAn + 1;
 589         if(ia <= 0) ia=1;
 590         if(ia > nofAnodes) ia = nofAnodes;
 591         aExpo     = (aStep*(ka-0.5)-aConst);
 592         if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
 593         else {
 594           Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
 595           anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
 596         }
 597         // index starts from 0
 598         index = iWing*nofAnodes+ia-1;
 599         if(anodeAmplitude){
 600           for(kt=jtmin; kt<=jtmax; kt++) {
 601             it = (kt-1)/nsplitTb+1;  // it starts from 1
 602             if(it<=0) it=1;
 603             if(it>fScaleSize*fMaxNofSamples)
 604               it = fScaleSize*fMaxNofSamples;
 605             tExpo    = (tStep*(kt-0.5)-tConst);
 606             if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
 607             else {
 608               Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
 609               timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep;
 610             }
 611             timeAmplitude *= nanoampToADC;
 612             //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
 613             Double_t charge = timeAmplitude;
 614             charge += fHitMap2->GetSignal(index,it-1);
 615             fHitMap2->SetHit(index, it-1, charge);
 616             fpList->AddSignal(index,it-1,itrack,ii-1,
 617                               mod->GetIndex(),timeAmplitude);
 618             fAnodeFire[index] = kTRUE;
 619           }  // end loop over time in window
 620         } // end if anodeAmplitude
 621       } // loop over anodes in window
 622     } // end loop over "sub-hits"
 623   } // end loop over hits
 624 }
 625
 626 //____________________________________________
 627 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) {
 628   // Adds a Digit.
 629   Int_t size = AliITSdigit::GetNTracks();
 630
 631   Int_t digits[3];
 632   Int_t * tracks = new Int_t[size];
 633   Int_t * hits = new Int_t[size];
 634   Float_t phys;
 635   Float_t * charges = new Float_t[size];
 636
 637   digits[0] = i;
 638   digits[1] = j;
 639   digits[2] = signalc;
 640
 641   AliITSpListItem *pItem = fpList->GetpListItem( i, j );
 642   if( pItem == 0 ) {
 643     phys = 0.0;
 644     for( Int_t l=0; l<size; l++ ) {
 645       tracks[l]  = 0;
 646       hits[l]    = 0;
 647       charges[l] = 0.0;
 648     }
 649   } else {
 650     Int_t idtrack =  pItem->GetTrack( 0 );
 651     if( idtrack >= 0 ) phys = pItem->GetSignal();
 652     else phys = 0.0;
 653
 654     for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
 655       tracks[l]  = pItem->GetTrack( l );
 656       hits[l]    = pItem->GetHit( l );
 657       charges[l] = pItem->GetSignal( l );
 658     }else{
 659       tracks[l]  = -3;
 660       hits[l]    = -1;
 661       charges[l] = 0.0;
 662     }// end for if
 663   }
 664
 665   fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale );
 666   delete [] tracks;
 667   delete [] hits;
 668   delete [] charges;
 669 }
 670 //______________________________________________________________________
 671 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) {
 672   // add baseline, noise, gain, electronics and ADC saturation effects
 673   // apply dead channels
 674
 675   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 676   Double_t baseline=0;
 677   Double_t noise=0;
 678   Double_t gain=0;
 679   Float_t contrib=0;
 680   Int_t i,k,kk;
 681   AliITSSimuParam* simpar = fDetType->GetSimuParam();
 682   Float_t maxadc = simpar->GetSDDMaxAdc();
 683   Int_t nGroup=fScaleSize;
 684   if(res->IsAMAt20MHz()){
 685     nGroup=fScaleSize/2;
 686   }
 687
 688   for (i=0;i<fNofMaps;i++) {
 689     if( !fAnodeFire[i] ) continue;
 690     baseline = res->GetBaseline(i);
 691     noise = res->GetNoise(i);
 692     gain = res->GetChannelGain(i)/fDetType->GetAverageGainSDD();
 693     if(res->IsBad()) gain=0.;
 694     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
 695     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 696       fInZR[k]  = fHitMap2->GetSignal(i,k);
 697       if(bAddGain) fInZR[k]*=gain;
 698       if( bAddNoise ) {
 699         contrib   = (baseline + noise*gRandom->Gaus());
 700         fInZR[k] += contrib;
 701       }
 702       fInZI[k]  = 0.;
 703     } // end for k
 704     if(!fDoFFT) {
 705       for(k=0; k<fMaxNofSamples; k++) {
 706         Double_t newcont = 0.;
 707         Double_t maxcont = 0.;
 708         for(kk=0;kk<fScaleSize;kk++) {
 709           newcont = fInZR[fScaleSize*k+kk];
 710           if(newcont > maxcont) maxcont = newcont;
 711         } // end for kk
 712         newcont = maxcont;
 713         if (newcont >= maxadc) newcont = maxadc -1;
 714         if(newcont >= baseline){
 715           Warning("","newcont=%f>=baseline=%f",newcont,baseline);
 716         } // end if
 717           // back to analog: ?
 718         fHitMap2->SetHit(i,k,newcont);
 719       }  // end for k
 720     }else{
 721       FastFourierTransform(&fInZR[0],&fInZI[0],1);
 722       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 723         Double_t rw = fElectronics->GetTraFunReal(k);
 724         Double_t iw = fElectronics->GetTraFunImag(k);
 725         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 726         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 727       } // end for k
 728       FastFourierTransform(&fOutZR[0],&fOutZI[0],-1);
 729       for(k=0; k<fMaxNofSamples; k++) {
 730         Double_t newcont1 = 0.;
 731         Double_t maxcont1 = 0.;
 732         for(kk=0;kk<nGroup;kk++) {
 733           newcont1 = fOutZR[fScaleSize*k+kk];
 734           if(newcont1 > maxcont1) maxcont1 = newcont1;
 735         } // end for kk
 736         newcont1 = maxcont1;
 737         if (newcont1 >= maxadc) newcont1 = maxadc -1;
 738         fHitMap2->SetHit(i,k,newcont1);
 739       } // end for k
 740     }
 741   } // end for i loop over anodes
 742   return;
 743 }
 744
 745 //______________________________________________________________________
 746 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 747     // function add the crosstalk effect to signal
 748     // temporal function, should be checked...!!!
 749
 750     // create and inizialice crosstalk map
 751     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 752     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 753     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 754     for( Int_t z=0; z<fNofMaps; z++ ) {
 755       Double_t baseline = calibr->GetBaseline(z);
 756         Bool_t on = kFALSE;
 757         Int_t tstart = 0;
 758         Int_t tstop = 0;
 759         Int_t nTsteps = 0;
 760
 761         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 762             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 763             if( fadc > baseline ) {
 764                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 765                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 766                     if( fadc1 < fadc ) continue;
 767                     on = kTRUE;
 768                     nTsteps = 0;
 769                     tstart = l;
 770                 }
 771                 nTsteps++;
 772             }
 773             else { // end fadc > baseline
 774                 if( on == kTRUE ) {
 775                     if( nTsteps > 2 ) {
 776                         tstop = l;
 777                         // make smooth derivative
 778                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 779                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 780                         for( Int_t i=tstart; i<tstop; i++ ) {
 781                             if( i > 2 && i < fMaxNofSamples-2 )
 782                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
 783                                     -0.1*fHitMap2->GetSignal( z,i-1 )
 784                                     +0.1*fHitMap2->GetSignal( z,i+1 )
 785                                     +0.2*fHitMap2->GetSignal( z,i+2 );
 786                         }
 787
 788                         // add crosstalk contribution to neibourg anodes
 789                         for( Int_t i=tstart; i<tstop; i++ ) {
 790                             Int_t anode = z - 1;
 791                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
 792                             Float_t ctktmp =  -dev[i1] * 0.25;
 793                             if( anode > 0 ) {
 794                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 795                             }
 796                             anode = z + 1;
 797                             if( anode < fNofMaps ) {
 798                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 799                             }
 800                         }
 801                         delete [] dev;
 802
 803                     } // if( nTsteps > 2 )
 804                     on = kFALSE;
 805                 }  // if( on == kTRUE )
 806             }  // else
 807         }
 808     }
 809
 810     for( Int_t a=0; a<fNofMaps; a++ )
 811         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
 812             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
 813             fHitMap2->SetHit( a, t, signal );
 814         }
 815
 816     delete [] ctk;
 817 }
 818
 819 //______________________________________________________________________
 820 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
 821     // To the 10 to 8 bit lossive compression.
 822     // code from Davide C. and Albert W.
 823
 824     if (signal < 128)  return signal;
 825     if (signal < 256)  return (128+((signal-128)>>1));
 826     if (signal < 512)  return (192+((signal-256)>>3));
 827     if (signal < 1024) return (224+((signal-512)>>4));
 828     return 0;
 829 }
 830 //______________________________________________________________________
 831 Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const {
 832   // Decompression from 8 to 10 bit
 833
 834   if (signal < 0 || signal > 255) {
 835     AliWarning(Form("Signal value %d out of range",signal));
 836     return 0;
 837   } // end if signal <0 || signal >255
 838
 839   if (signal < 128) return signal;
 840   if (signal < 192) {
 841     if (TMath::Odd(signal)) return (128+((signal-128)<<1));
 842     else  return (128+((signal-128)<<1)+1);
 843   } // end if signal < 192
 844   if (signal < 224) {
 845     if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3);
 846     else  return (256+((signal-192)<<3)+4);
 847   } // end if signal < 224
 848   if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7);
 849   return (512+((signal-224)<<4)+8);
 850 }
 851 //______________________________________________________________________
 852 void AliITSsimulationSDD::Compress2D(){
 853   // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10
 854   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 855   for (Int_t iWing=0; iWing<2; iWing++) {
 856     Int_t tL=res->GetZSLowThreshold(iWing);
 857     Int_t tH=res->GetZSHighThreshold(iWing);
 858     for (Int_t i=0; i<fNofMaps/2; i++) {
 859       Int_t ian=i+iWing*fNofMaps/2;
 860       if( !fAnodeFire[ian] ) continue;
 861       for (Int_t itb=0; itb<fMaxNofSamples; itb++) {
 862         Int_t nLow=0, nHigh=0;
 863         Float_t cC=fHitMap2->GetSignal(ian,itb);
 864         if(cC<=tL) continue;
 865         nLow++; // cC is greater than tL
 866         if(cC>tH) nHigh++;
 867         //                     N
 868         // Get "quintuple":   WCE
 869         //                     S
 870         Float_t wW=0.;
 871         if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1);
 872         if(wW>tL) nLow++;
 873         if(wW>tH) nHigh++;
 874         Float_t eE=0.;
 875         if(itb<fMaxNofSamples-1) eE=fHitMap2->GetSignal(ian,itb+1);
 876         if(eE>tL) nLow++;
 877         if(eE>tH) nHigh++;
 878         Float_t nN=0.;
 879         if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb);
 880         if(nN>tL) nLow++;
 881         if(nN>tH) nHigh++;
 882         Float_t sS=0.;
 883         if(i>0) sS=fHitMap2->GetSignal(ian-1,itb);
 884         if(sS>tL) nLow++;
 885         if(sS>tH) nHigh++;
 886
 887         if(nLow>=2 && nHigh>=1){
 888           Int_t signal=(Int_t)cC;
 889           Int_t signalc = Convert10to8(signal);
 890           Int_t signale = Convert8to10(signalc);
 891           signalc-=tL; // subtract low threshold after 10 to 8 bit compression
 892           if(signalc>=4) AddDigit(ian,itb,signalc,signale);  // store C
 893         }
 894       }
 895     }
 896   }
 897 }
 898
 899
 900 //______________________________________________________________________
 901 void AliITSsimulationSDD::StoreAllDigits(){
 902   // store digits for non-zero-suppressed data
 903   for (Int_t ian=0; ian<fNofMaps; ian++) {
 904     for (Int_t itb=0; itb<fMaxNofSamples; itb++){
 905       Int_t signal=(Int_t)(fHitMap2->GetSignal(ian,itb));
 906       Int_t signalc = Convert10to8(signal);
 907       Int_t signale = Convert8to10(signalc);
 908       AddDigit(ian,itb,signalc,signale);
 909     }
 910   }
 911 }
 912 //______________________________________________________________________
 913 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
 914   // Creates histograms of maps for debugging
 915   Int_t i;
 916
 917   fHis=new TObjArray(fNofMaps);
 918   for (i=0;i<fNofMaps;i++) {
 919     TString sddName;
 920     sddName.Form("sdd_%d",i+1);
 921     fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
 922                          0.,(Float_t) scale*fMaxNofSamples), i);
 923   } // end for i
 924 }
 925 //______________________________________________________________________
 926 void AliITSsimulationSDD::FillHistograms(){
 927     // fill 1D histograms from map
 928
 929     if (!fHis) return;
 930
 931     for( Int_t i=0; i<fNofMaps; i++) {
 932         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
 933         Int_t nsamples = hist->GetNbinsX();
 934         for( Int_t j=0; j<nsamples; j++) {
 935             Double_t signal=fHitMap2->GetSignal(i,j);
 936             hist->Fill((Float_t)j,signal);
 937         } // end for j
 938     } // end for i
 939 }
 940 //______________________________________________________________________
 941 void AliITSsimulationSDD::ResetHistograms(){
 942     // Reset histograms for this detector
 943     Int_t i;
 944
 945     for (i=0;i<fNofMaps;i++ ) {
 946         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
 947     } // end for i
 948 }
 949 //______________________________________________________________________
 950 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
 951     // Fills a histogram from a give anode.
 952
 953     if (!fHis) return 0;
 954
 955     if(wing <=0 || wing > 2) {
 956         Warning("GetAnode","Wrong wing number: %d",wing);
 957         return NULL;
 958     } // end if wing <=0 || wing >2
 959     if(anode <=0 || anode > fNofMaps/2) {
 960         Warning("GetAnode","Wrong anode number: %d",anode);
 961         return NULL;
 962     } // end if ampde <=0 || andoe > fNofMaps/2
 963
 964     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
 965     return (TH1F*)(fHis->At(index));
 966 }
 967 //______________________________________________________________________
 968 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
 969     // Writes the histograms to a file
 970
 971     if (!fHis) return;
 972
 973     hfile->cd();
 974     Int_t i;
 975     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
 976     return;
 977 }
 978 //______________________________________________________________________
 979 void AliITSsimulationSDD::WriteSDigits(){
 980     // Fills the Summable digits Tree
 981     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 982
 983     for( Int_t i=0; i<fNofMaps; i++ ) {
 984         if( !fAnodeFire[i] ) continue;
 985         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 986             Double_t sig = fHitMap2->GetSignal( i, j );
 987             if( sig > 0.2 ) {
 988                 Int_t jdx = j*fScaleSize;
 989                 Int_t index = fpList->GetHitIndex( i, j );
 990                 AliITSpListItem pItemTmp2( fModule, index, 0. );
 991                 // put the fScaleSize analog digits in only one
 992                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 993                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
 994                     if( pItemTmp == 0 ) continue;
 995                     pItemTmp2.Add( pItemTmp );
 996                 }
 997                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
 998                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
 999                 aliITS->AddSumDigit( pItemTmp2 );
1000             } // end if (sig > 0.2)
1001         }
1002     }
1003     return;
1004 }
1005 //______________________________________________________________________
1006 void AliITSsimulationSDD::PrintStatus() const {
1007     // Print SDD simulation Parameters
1008
1009     cout << "**************************************************" << endl;
1010     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1011     cout << "**************************************************" << endl;
1012     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1013     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1014     cout << "Number of Anodes used: " << fNofMaps << endl;
1015     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1016     cout << "Scale size factor: " << fScaleSize << endl;
1017     cout << "**************************************************" << endl;
1018 }