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