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