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