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