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