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