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