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 "AliITSpList.h"
  36 #include "AliITSCalibrationSDD.h"
  37 #include "AliITSsegmentationSDD.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     SetScaleFourier();
  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     SetScaleFourier();
 148     SetPerpendTracksFlag();
 149     SetCrosstalkFlag();
 150     SetDoFFT();
 151
 152     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 153     if(seg->Npx()==128) fScaleSize=8;
 154     AliITSSimuParam* simpar = fDetType->GetSimuParam();
 155     fpList = new AliITSpList( seg->Npz(),
 156                               fScaleSize*seg->Npx() );
 157     fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
 158     fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
 159     fHitMap2 = fHitSigMap2;
 160
 161     fNofMaps = seg->Npz();
 162     fMaxNofSamples = seg->Npx();
 163     fAnodeFire = new Bool_t [fNofMaps];
 164
 165     Float_t sddWidth  = seg->Dz();
 166     Float_t anodePitch = seg->Dpz(0);
 167     Double_t timeStep  = (Double_t)seg->Dpx(0);
 168
 169     if(anodePitch*(fNofMaps/2) > sddWidth) {
 170         Warning("AliITSsimulationSDD",
 171                 "Too many anodes %d or too big pitch %f \n",
 172                 fNofMaps/2,anodePitch);
 173     } // end if
 174
 175
 176     fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
 177                                     simpar->GetSDDElectronics());
 178
 179
 180     fITS       = (AliITS*)gAlice->GetModule("ITS");
 181
 182     fInZR  = new Double_t [fScaleSize*fMaxNofSamples];
 183     fInZI  = new Double_t [fScaleSize*fMaxNofSamples];
 184     fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
 185     fOutZI = new Double_t [fScaleSize*fMaxNofSamples];
 186 }
 187 //______________________________________________________________________
 188 AliITSsimulationSDD::~AliITSsimulationSDD() {
 189     // destructor
 190
 191     //    delete fpList;
 192     delete fHitSigMap2;
 193     delete fHitNoiMap2;
 194     delete fElectronics;
 195
 196     fITS = 0;
 197
 198     if (fHis) {
 199         fHis->Delete();
 200         delete fHis;
 201     } // end if fHis
 202     if(fInZR)  delete [] fInZR;
 203     if(fInZI)  delete [] fInZI;
 204     if(fOutZR) delete [] fOutZR;
 205     if(fOutZI) delete [] fOutZI;
 206     if(fAnodeFire) delete [] fAnodeFire;
 207 }
 208 //______________________________________________________________________
 209 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
 210     // create maps to build the lists of tracks for each summable digit
 211     fModule = module;
 212     fEvent  = event;
 213     ClearMaps();
 214     memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);
 215 }
 216 //______________________________________________________________________
 217 void AliITSsimulationSDD::ClearMaps() {
 218     // clear maps
 219     fpList->ClearMap();
 220     fHitSigMap2->ClearMap();
 221     fHitNoiMap2->ClearMap();
 222 }
 223 //______________________________________________________________________
 224 void AliITSsimulationSDD::FastFourierTransform(Double_t *real,
 225                           Double_t *imag,Int_t direction) {
 226     // Do a Fast Fourier Transform
 227
 228     Int_t samples = fElectronics->GetSamples();
 229     Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
 230     Int_t m1 = samples;
 231     Int_t m  = samples/2;
 232     Int_t m2 = samples/m1;
 233     Int_t i,j,k;
 234     for(i=1; i<=l; i++) {
 235         for(j=0; j<samples; j += m1) {
 236             Int_t p = 0;
 237             for(k=j; k<= j+m-1; k++) {
 238                 Double_t wsr = fElectronics->GetWeightReal(p);
 239                 Double_t wsi = fElectronics->GetWeightImag(p);
 240                 if(direction == -1) wsi = -wsi;
 241                 Double_t xr = *(real+k+m);
 242                 Double_t xi = *(imag+k+m);
 243                 *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi);
 244                 *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr);
 245                 *(real+k) += xr;
 246                 *(imag+k) += xi;
 247                 p += m2;
 248             } // end for k
 249         } // end for j
 250         m1 = m;
 251         m /= 2;
 252         m2 += m2;
 253     } // end for i
 254     for(j=0; j<samples; j++) {
 255         Int_t j1 = j;
 256         Int_t p = 0;
 257         Int_t i1;
 258         for(i1=1; i1<=l; i1++) {
 259             Int_t j2 = j1;
 260             j1 /= 2;
 261             p = p + p + j2 - j1 - j1;
 262         } // end for i1
 263         if(p >= j) {
 264             Double_t xr = *(real+j);
 265             Double_t xi = *(imag+j);
 266             *(real+j) = *(real+p);
 267             *(imag+j) = *(imag+p);
 268             *(real+p) = xr;
 269             *(imag+p) = xi;
 270         } // end if p>=j
 271     } // end for j
 272     if(direction == -1) {
 273         for(i=0; i<samples; i++) {
 274             *(real+i) /= samples;
 275             *(imag+i) /= samples;
 276         } // end for i
 277     } // end if direction == -1
 278     return;
 279 }
 280
 281 //______________________________________________________________________
 282 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 283     // digitize module using the "slow" detector simulator creating
 284     // summable digits.
 285
 286     TObjArray *fHits = mod->GetHits();
 287     Int_t nhits      = fHits->GetEntriesFast();
 288     if( !nhits ) return;
 289
 290     InitSimulationModule( md, ev );
 291     HitsToAnalogDigits( mod );  // fills fHitMap2 which is = fHitSigmap2
 292     ChargeToSignal( fModule,kFALSE,kTRUE ); // - Process signal adding gain without adding noise
 293     fHitMap2 = fHitNoiMap2;   // - Swap to noise map
 294     ChargeToSignal( fModule,kTRUE,kFALSE );  // - Process only noise
 295     fHitMap2 = fHitSigMap2;   // - Return to signal map
 296     WriteSDigits();
 297     ClearMaps();
 298 }
 299 //______________________________________________________________________
 300 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
 301                                                Int_t mask ) {
 302     // Add Summable digits to module maps.
 303     AliITSSimuParam* simpar = fDetType->GetSimuParam();
 304     Int_t    nItems = pItemArray->GetEntries();
 305     Double_t maxadc = simpar->GetSDDMaxAdc();
 306     Bool_t sig = kFALSE;
 307
 308     // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
 309     for( Int_t i=0; i<nItems; i++ ) {
 310         AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
 311         if( pItem->GetModule() != fModule ) {
 312             Error( "AliITSsimulationSDD","Error reading, SDigits module "
 313                    "%d != current module %d: exit",
 314                    pItem->GetModule(), fModule );
 315             return sig;
 316         } // end if
 317
 318         if(pItem->GetSignal()>0.0 ) sig = kTRUE;
 319
 320         fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
 321         AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
 322         Double_t sigAE = pItem2->GetSignalAfterElect();
 323         if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
 324         Int_t ia;
 325         Int_t it;
 326         fpList->GetMapIndex( pItem->GetIndex(), ia, it );
 327         fHitMap2->SetHit( ia, it, sigAE );
 328         fAnodeFire[ia] = kTRUE;
 329     }
 330     return sig;
 331 }
 332 //______________________________________________________________________
 333 void AliITSsimulationSDD::FinishSDigitiseModule() {
 334     // digitize module using the "slow" detector simulator from
 335     // the sum of summable digits.
 336     FinishDigits() ;
 337     ClearMaps();
 338 }
 339 //______________________________________________________________________
 340 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 341     // create maps to build the lists of tracks for each digit
 342
 343     TObjArray *fHits = mod->GetHits();
 344     Int_t nhits      = fHits->GetEntriesFast();
 345
 346     InitSimulationModule( md, ev );
 347     if( !nhits ) return;
 348
 349     HitsToAnalogDigits( mod );
 350     ChargeToSignal( fModule,kTRUE,kTRUE );  // process signal + noise
 351
 352     for( Int_t i=0; i<fNofMaps; i++ ) {
 353         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 354             Int_t jdx = j*fScaleSize;
 355             Int_t index = fpList->GetHitIndex( i, j );
 356             AliITSpListItem pItemTmp2( fModule, index, 0. );
 357             // put the fScaleSize analog digits in only one
 358             for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 359                 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
 360                 if( pItemTmp == 0 ) continue;
 361                 pItemTmp2.Add( pItemTmp );
 362             }
 363             fpList->DeleteHit( i, j );
 364             fpList->AddItemTo( 0, &pItemTmp2 );
 365         }
 366     }
 367     FinishDigits();
 368     ClearMaps();
 369 }
 370 //______________________________________________________________________
 371 void AliITSsimulationSDD::FinishDigits() {
 372     // introduce the electronics effects and do zero-suppression if required
 373
 374     if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
 375
 376     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 377     Bool_t isZeroSupp = res->GetZeroSupp();
 378     if (isZeroSupp) Compress2D();
 379     else StoreAllDigits();
 380 }
 381 //______________________________________________________________________
 382 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
 383     // create maps to build the lists of tracks for each digit
 384   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 385   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 386   AliITSSimuParam* simpar = fDetType->GetSimuParam();
 387   TObjArray *hits     = mod->GetHits();
 388     Int_t      nhits    = hits->GetEntriesFast();
 389
 390     //    Int_t      arg[6]   = {0,0,0,0,0,0};
 391     Int_t     nofAnodes  = fNofMaps/2;
 392     Double_t  sddLength  = seg->Dx();
 393     Double_t  sddWidth   = seg->Dz();
 394     Double_t  anodePitch = seg->Dpz(0);
 395     Double_t  timeStep   = seg->Dpx(0);
 396     Double_t  driftSpeed ;  // drift velocity (anode dependent)
 397     Double_t  nanoampToADC       = simpar->GetSDDMaxAdc()/simpar->GetSDDDynamicRange(); //   maxadc/topValue;
 398     Double_t  cHloss     = simpar->GetSDDChargeLoss();
 399     Float_t   dfCoeff, s1;
 400     simpar->GetSDDDiffCoeff(dfCoeff,s1); // Signal 2d Shape
 401     Double_t  eVpairs    = simpar->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
 402     Double_t  nsigma     = simpar->GetNSigmaIntegration(); //
 403     Int_t     nlookups   = simpar->GetGausNLookUp();       //
 404     Float_t   jitter     = simpar->GetSDDJitterError(); //
 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\n");
 461       }
 462       depEnergy  *= kconv;
 463
 464       // scale path to simulate a perpendicular track
 465       // continue if the particle did not lose energy
 466       // passing through detector
 467       if (!depEnergy) {
 468         AliDebug(1,
 469                  Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
 470                       itrack,ii,mod->GetIndex()));
 471         continue;
 472       } // end if !depEnergy
 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         timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1); // time bin in range 1-256 !!!
 516         if(timeSample > fScaleSize*fMaxNofSamples) {
 517           AliWarning(Form("Wrong Time Sample: %e",timeSample));
 518           continue;
 519         } // end if timeSample > fScaleSize*fMaxNofSamples
 520         if(zAnode>nofAnodes) zAnode-=nofAnodes;  // to have the anode number between 0. and 256.
 521         if(zAnode*anodePitch > sddWidth || zAnode*anodePitch < 0.)
 522           AliWarning(Form("Exceeding sddWidth=%e Z = %e",sddWidth,zAnode*anodePitch));
 523         iAnode = (Int_t) (1.+zAnode); // iAnode in range 1-256 !!!!
 524         if(iAnode < 1 || iAnode > nofAnodes) {
 525           AliWarning(Form("Wrong iAnode: 1<%d>%d  (xanode=%e)",iAnode,nofAnodes, zAnode));
 526           continue;
 527         } // end if iAnode < 1 || iAnode > nofAnodes
 528
 529         // store straight away the particle position in the array
 530         // of particles and take idhit=ii only when part is entering (this
 531         // requires FillModules() in the macro for analysis) :
 532
 533         // Sigma along the anodes for track segment.
 534         sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
 535         sigT       = sigA/driftSpeed;
 536         // Peak amplitude in nanoAmpere
 537         amplitude  = fScaleSize*160.*depEnergy/
 538           (timeStep*eVpairs*2.*acos(-1.));
 539         chargeloss = 1.-cHloss*driftPath/1000.;
 540         amplitude *= chargeloss;
 541         width  = 2.*nsigma/(nlookups-1);
 542         // Spread the charge
 543         nsplitAn = 4;
 544         nsplitTb=4;
 545         aStep  = anodePitch/(nsplitAn*sigA);
 546         aConst = zAnode*anodePitch/sigA;
 547         tStep  = timeStep/(nsplitTb*fScaleSize*sigT);
 548         tConst = drTime/sigT;
 549         // Define SDD window corresponding to the hit
 550         anodeWindow = (Int_t)(nsigma*sigA/anodePitch+1);
 551         timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
 552         jamin = (iAnode - anodeWindow - 2)*nsplitAn+1;
 553         jamax = (iAnode + anodeWindow + 2)*nsplitAn;
 554         if(jamin <= 0) jamin = 1;
 555         if(jamax > nofAnodes*nsplitAn)
 556           jamax = nofAnodes*nsplitAn;
 557         // jtmin and jtmax are Hard-wired
 558         jtmin = (Int_t)(timeSample-timeWindow-2)*nsplitTb+1;
 559         jtmax = (Int_t)(timeSample+timeWindow+2)*nsplitTb;
 560         if(jtmin <= 0) jtmin = 1;
 561         if(jtmax > fScaleSize*fMaxNofSamples*nsplitTb)
 562           jtmax = fScaleSize*fMaxNofSamples*nsplitTb;
 563         // Spread the charge in the anode-time window
 564         for(ka=jamin; ka <=jamax; ka++) {
 565           ia = (ka-1)/nsplitAn + 1;
 566           if(ia <= 0) ia=1;
 567           if(ia > nofAnodes) ia = nofAnodes;
 568           aExpo     = (aStep*(ka-0.5)-aConst);
 569           if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
 570           else {
 571             Int_t theBin = (Int_t) ((aExpo+nsigma)/width+0.5);
 572             anodeAmplitude = amplitude*simpar->GetGausLookUp(theBin);
 573           }
 574           // index starts from 0
 575           index = iWing*nofAnodes+ia-1;
 576           if(anodeAmplitude){
 577             for(kt=jtmin; kt<=jtmax; kt++) {
 578               it = (kt-1)/nsplitTb+1;  // it starts from 1
 579               if(it<=0) it=1;
 580               if(it>fScaleSize*fMaxNofSamples)
 581                 it = fScaleSize*fMaxNofSamples;
 582               tExpo    = (tStep*(kt-0.5)-tConst);
 583               if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
 584               else {
 585                 Int_t theBin = (Int_t) ((tExpo+nsigma)/width+0.5);
 586                 timeAmplitude = anodeAmplitude*simpar->GetGausLookUp(theBin)*aStep*tStep;
 587               }
 588               timeAmplitude *= nanoampToADC;
 589               //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
 590               Double_t charge = timeAmplitude;
 591               charge += fHitMap2->GetSignal(index,it-1);
 592               fHitMap2->SetHit(index, it-1, charge);
 593               fpList->AddSignal(index,it-1,itrack,ii-1,
 594                                 mod->GetIndex(),timeAmplitude);
 595               fAnodeFire[index] = kTRUE;
 596             }  // end loop over time in window
 597           } // end if anodeAmplitude
 598         } // loop over anodes in window
 599       } // end loop over "sub-hits"
 600     } // end loop over hits
 601 }
 602
 603 //____________________________________________
 604 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signalc, Int_t signale) {
 605   // Adds a Digit.
 606   Int_t size = AliITSdigit::GetNTracks();
 607
 608   Int_t digits[3];
 609   Int_t * tracks = new Int_t[size];
 610   Int_t * hits = new Int_t[size];
 611   Float_t phys;
 612   Float_t * charges = new Float_t[size];
 613
 614   digits[0] = i;
 615   digits[1] = j;
 616   digits[2] = signalc;
 617
 618   AliITSpListItem *pItem = fpList->GetpListItem( i, j );
 619   if( pItem == 0 ) {
 620     phys = 0.0;
 621     for( Int_t l=0; l<size; l++ ) {
 622       tracks[l]  = 0;
 623       hits[l]    = 0;
 624       charges[l] = 0.0;
 625     }
 626   } else {
 627     Int_t idtrack =  pItem->GetTrack( 0 );
 628     if( idtrack >= 0 ) phys = pItem->GetSignal();
 629     else phys = 0.0;
 630
 631     for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
 632       tracks[l]  = pItem->GetTrack( l );
 633       hits[l]    = pItem->GetHit( l );
 634       charges[l] = pItem->GetSignal( l );
 635     }else{
 636       tracks[l]  = -3;
 637       hits[l]    = -1;
 638       charges[l] = 0.0;
 639     }// end for if
 640   }
 641
 642   fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges, signale );
 643   delete [] tracks;
 644   delete [] hits;
 645   delete [] charges;
 646 }
 647 //______________________________________________________________________
 648 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise, Bool_t bAddGain) {
 649   // add baseline, noise, gain, electronics and ADC saturation effects
 650   // apply dead channels
 651
 652   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 653   Double_t baseline=0;
 654   Double_t noise=0;
 655   Double_t gain=0;
 656   Float_t contrib=0;
 657   Int_t i,k,kk;
 658   AliITSSimuParam* simpar = fDetType->GetSimuParam();
 659   Float_t maxadc = simpar->GetSDDMaxAdc();
 660   Int_t nGroup=fScaleSize;
 661   if(res->IsAMAt20MHz()){
 662     nGroup=fScaleSize/2;
 663   }
 664
 665   for (i=0;i<fNofMaps;i++) {
 666     if( !fAnodeFire[i] ) continue;
 667     baseline = res->GetBaseline(i);
 668     noise = res->GetNoise(i);
 669     gain = res->GetChannelGain(i);
 670     if(res->IsBad()) gain=0.;
 671     if( res->IsChipBad(res->GetChip(i)) )gain=0.;
 672     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 673       fInZR[k]  = fHitMap2->GetSignal(i,k);
 674       if(bAddGain) fInZR[k]*=gain;
 675       if( bAddNoise ) {
 676         contrib   = (baseline + noise*gRandom->Gaus());
 677         fInZR[k] += contrib;
 678       }
 679       fInZI[k]  = 0.;
 680     } // end for k
 681     if(!fDoFFT) {
 682       for(k=0; k<fMaxNofSamples; k++) {
 683         Double_t newcont = 0.;
 684         Double_t maxcont = 0.;
 685         for(kk=0;kk<fScaleSize;kk++) {
 686           newcont = fInZR[fScaleSize*k+kk];
 687           if(newcont > maxcont) maxcont = newcont;
 688         } // end for kk
 689         newcont = maxcont;
 690         if (newcont >= maxadc) newcont = maxadc -1;
 691         if(newcont >= baseline){
 692           Warning("","newcont=%d>=baseline=%d",newcont,baseline);
 693         } // end if
 694           // back to analog: ?
 695         fHitMap2->SetHit(i,k,newcont);
 696       }  // end for k
 697     }else{
 698       FastFourierTransform(&fInZR[0],&fInZI[0],1);
 699       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 700         Double_t rw = fElectronics->GetTraFunReal(k);
 701         Double_t iw = fElectronics->GetTraFunImag(k);
 702         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 703         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 704       } // end for k
 705       FastFourierTransform(&fOutZR[0],&fOutZI[0],-1);
 706       for(k=0; k<fMaxNofSamples; k++) {
 707         Double_t newcont1 = 0.;
 708         Double_t maxcont1 = 0.;
 709         for(kk=0;kk<nGroup;kk++) {
 710           newcont1 = fOutZR[fScaleSize*k+kk];
 711           if(newcont1 > maxcont1) maxcont1 = newcont1;
 712         } // end for kk
 713         newcont1 = maxcont1;
 714         if (newcont1 >= maxadc) newcont1 = maxadc -1;
 715         fHitMap2->SetHit(i,k,newcont1);
 716       } // end for k
 717     }
 718   } // end for i loop over anodes
 719   return;
 720 }
 721
 722 //______________________________________________________________________
 723 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 724     // function add the crosstalk effect to signal
 725     // temporal function, should be checked...!!!
 726
 727     // create and inizialice crosstalk map
 728     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 729     if( ctk == NULL ) {
 730         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
 731         return;
 732     }
 733     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 734     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 735     for( Int_t z=0; z<fNofMaps; z++ ) {
 736       Double_t baseline = calibr->GetBaseline(z);
 737         Bool_t on = kFALSE;
 738         Int_t tstart = 0;
 739         Int_t tstop = 0;
 740         Int_t nTsteps = 0;
 741
 742         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 743             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 744             if( fadc > baseline ) {
 745                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 746                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 747                     if( fadc1 < fadc ) continue;
 748                     on = kTRUE;
 749                     nTsteps = 0;
 750                     tstart = l;
 751                 }
 752                 nTsteps++;
 753             }
 754             else { // end fadc > baseline
 755                 if( on == kTRUE ) {
 756                     if( nTsteps > 2 ) {
 757                         tstop = l;
 758                         // make smooth derivative
 759                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 760                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 761                         if( ctk == NULL ) {
 762                             Error( "ApplyCrosstalk",
 763                                    "no memory for temporal array: exit \n" );
 764                             return;
 765                         }
 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>=3 && 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           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("sdd_");
 906         Char_t candNum[4];
 907         sprintf(candNum,"%d",i+1);
 908         sddName.Append(candNum);
 909         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
 910                              0.,(Float_t) scale*fMaxNofSamples), i);
 911     } // end for i
 912 }
 913 //______________________________________________________________________
 914 void AliITSsimulationSDD::FillHistograms(){
 915     // fill 1D histograms from map
 916
 917     if (!fHis) return;
 918
 919     for( Int_t i=0; i<fNofMaps; i++) {
 920         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
 921         Int_t nsamples = hist->GetNbinsX();
 922         for( Int_t j=0; j<nsamples; j++) {
 923             Double_t signal=fHitMap2->GetSignal(i,j);
 924             hist->Fill((Float_t)j,signal);
 925         } // end for j
 926     } // end for i
 927 }
 928 //______________________________________________________________________
 929 void AliITSsimulationSDD::ResetHistograms(){
 930     // Reset histograms for this detector
 931     Int_t i;
 932
 933     for (i=0;i<fNofMaps;i++ ) {
 934         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
 935     } // end for i
 936 }
 937 //______________________________________________________________________
 938 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
 939     // Fills a histogram from a give anode.
 940
 941     if (!fHis) return 0;
 942
 943     if(wing <=0 || wing > 2) {
 944         Warning("GetAnode","Wrong wing number: %d",wing);
 945         return NULL;
 946     } // end if wing <=0 || wing >2
 947     if(anode <=0 || anode > fNofMaps/2) {
 948         Warning("GetAnode","Wrong anode number: %d",anode);
 949         return NULL;
 950     } // end if ampde <=0 || andoe > fNofMaps/2
 951
 952     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
 953     return (TH1F*)(fHis->At(index));
 954 }
 955 //______________________________________________________________________
 956 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
 957     // Writes the histograms to a file
 958
 959     if (!fHis) return;
 960
 961     hfile->cd();
 962     Int_t i;
 963     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
 964     return;
 965 }
 966 //______________________________________________________________________
 967 void AliITSsimulationSDD::WriteSDigits(){
 968     // Fills the Summable digits Tree
 969     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 970
 971     for( Int_t i=0; i<fNofMaps; i++ ) {
 972         if( !fAnodeFire[i] ) continue;
 973         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 974             Double_t sig = fHitMap2->GetSignal( i, j );
 975             if( sig > 0.2 ) {
 976                 Int_t jdx = j*fScaleSize;
 977                 Int_t index = fpList->GetHitIndex( i, j );
 978                 AliITSpListItem pItemTmp2( fModule, index, 0. );
 979                 // put the fScaleSize analog digits in only one
 980                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 981                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
 982                     if( pItemTmp == 0 ) continue;
 983                     pItemTmp2.Add( pItemTmp );
 984                 }
 985                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
 986                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
 987                 aliITS->AddSumDigit( pItemTmp2 );
 988             } // end if (sig > 0.2)
 989         }
 990     }
 991     return;
 992 }
 993 //______________________________________________________________________
 994 void AliITSsimulationSDD::PrintStatus() const {
 995     // Print SDD simulation Parameters
 996
 997     cout << "**************************************************" << endl;
 998     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
 999     cout << "**************************************************" << endl;
1000     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1001     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1002     cout << "Number of Anodes used: " << fNofMaps << endl;
1003     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1004     cout << "Scale size factor: " << fScaleSize << endl;
1005     cout << "**************************************************" << endl;
1006 }