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