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