ITS/AliITSsimulationSDD.cxx

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