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