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 signalc, Int_t signale) {
 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] = signalc;
 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, signale );
 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
 753     // create and inizialice crosstalk map
 754     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 755     if( ctk == NULL ) {
 756         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
 757         return;
 758     }
 759     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 760     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 761     for( Int_t z=0; z<fNofMaps; z++ ) {
 762       Double_t baseline = calibr->GetBaseline(z);
 763         Bool_t on = kFALSE;
 764         Int_t tstart = 0;
 765         Int_t tstop = 0;
 766         Int_t nTsteps = 0;
 767
 768         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 769             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 770             if( fadc > baseline ) {
 771                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 772                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 773                     if( fadc1 < fadc ) continue;
 774                     on = kTRUE;
 775                     nTsteps = 0;
 776                     tstart = l;
 777                 }
 778                 nTsteps++;
 779             }
 780             else { // end fadc > baseline
 781                 if( on == kTRUE ) {
 782                     if( nTsteps > 2 ) {
 783                         tstop = l;
 784                         // make smooth derivative
 785                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 786                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 787                         if( ctk == NULL ) {
 788                             Error( "ApplyCrosstalk",
 789                                    "no memory for temporal array: exit \n" );
 790                             return;
 791                         }
 792                         for( Int_t i=tstart; i<tstop; i++ ) {
 793                             if( i > 2 && i < fMaxNofSamples-2 )
 794                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
 795                                     -0.1*fHitMap2->GetSignal( z,i-1 )
 796                                     +0.1*fHitMap2->GetSignal( z,i+1 )
 797                                     +0.2*fHitMap2->GetSignal( z,i+2 );
 798                         }
 799
 800                         // add crosstalk contribution to neibourg anodes
 801                         for( Int_t i=tstart; i<tstop; i++ ) {
 802                             Int_t anode = z - 1;
 803                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
 804                             Float_t ctktmp =  -dev[i1] * 0.25;
 805                             if( anode > 0 ) {
 806                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 807                             }
 808                             anode = z + 1;
 809                             if( anode < fNofMaps ) {
 810                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 811                             }
 812                         }
 813                         delete [] dev;
 814
 815                     } // if( nTsteps > 2 )
 816                     on = kFALSE;
 817                 }  // if( on == kTRUE )
 818             }  // else
 819         }
 820     }
 821
 822     for( Int_t a=0; a<fNofMaps; a++ )
 823         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
 824             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
 825             fHitMap2->SetHit( a, t, signal );
 826         }
 827
 828     delete [] ctk;
 829 }
 830
 831 //______________________________________________________________________
 832 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
 833     // To the 10 to 8 bit lossive compression.
 834     // code from Davide C. and Albert W.
 835
 836     if (signal < 128)  return signal;
 837     if (signal < 256)  return (128+((signal-128)>>1));
 838     if (signal < 512)  return (192+((signal-256)>>3));
 839     if (signal < 1024) return (224+((signal-512)>>4));
 840     return 0;
 841 }
 842 //______________________________________________________________________
 843 Int_t AliITSsimulationSDD::Convert8to10(Int_t signal) const {
 844   // Decompression from 8 to 10 bit
 845
 846   if (signal < 0 || signal > 255) {
 847     AliWarning(Form("Signal value %d out of range",signal));
 848     return 0;
 849   } // end if signal <0 || signal >255
 850
 851   if (signal < 128) return signal;
 852   if (signal < 192) {
 853     if (TMath::Odd(signal)) return (128+((signal-128)<<1));
 854     else  return (128+((signal-128)<<1)+1);
 855   } // end if signal < 192
 856   if (signal < 224) {
 857     if (TMath::Odd(signal)) return (256+((signal-192)<<3)+3);
 858     else  return (256+((signal-192)<<3)+4);
 859   } // end if signal < 224
 860   if (TMath::Odd(signal)) return (512+((signal-224)<<4)+7);
 861   return (512+((signal-224)<<4)+8);
 862 }
 863 //______________________________________________________________________
 864 void AliITSsimulationSDD::Compress2D(){
 865   // 2D zero-suppression algorithm as described in ALICE-INT-1999-28 V10
 866   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 867   for (Int_t iWing=0; iWing<2; iWing++) {
 868     Int_t tL=res->GetZSLowThreshold(iWing);
 869     Int_t tH=res->GetZSHighThreshold(iWing);
 870     for (Int_t i=0; i<fNofMaps/2; i++) {
 871       Int_t ian=i+iWing*fNofMaps/2;
 872       if( !fAnodeFire[ian] ) continue;
 873       for (Int_t itb=0; itb<fMaxNofSamples; itb++) {
 874         Int_t nLow=0, nHigh=0;
 875         Float_t cC=fHitMap2->GetSignal(ian,itb);
 876         if(cC<=tL) continue;
 877         nLow++; // cC is greater than tL
 878         if(cC>tH) nHigh++;
 879         //                     N
 880         // Get "quintuple":   WCE
 881         //                     S
 882         Float_t wW=0.;
 883         if(itb>0) wW=fHitMap2->GetSignal(ian,itb-1);
 884         if(wW>tL) nLow++;
 885         if(wW>tH) nHigh++;
 886         Float_t eE=0.;
 887         if(itb<fMaxNofSamples-1) eE=fHitMap2->GetSignal(ian,itb+1);
 888         if(eE>tL) nLow++;
 889         if(eE>tH) nHigh++;
 890         Float_t nN=0.;
 891         if(i<(fNofMaps/2-1)) nN=fHitMap2->GetSignal(ian+1,itb);
 892         if(nN>tL) nLow++;
 893         if(nN>tH) nHigh++;
 894         Float_t sS=0.;
 895         if(i>0) sS=fHitMap2->GetSignal(ian-1,itb);
 896         if(sS>tL) nLow++;
 897         if(sS>tH) nHigh++;
 898
 899         if(nLow>=3 && nHigh>=1){
 900           Int_t signal=(Int_t)cC;
 901           Int_t signalc = Convert10to8(signal);
 902           Int_t signale = Convert8to10(signalc);
 903           signalc-=tL; // subtract low threshold after 10 to 8 bit compression
 904           AddDigit(ian,itb,signalc,signale);  // store C
 905         }
 906       }
 907     }
 908   }
 909 }
 910
 911
 912 //______________________________________________________________________
 913 void AliITSsimulationSDD::StoreAllDigits(){
 914   // store digits for non-zero-suppressed data
 915   for (Int_t ian=0; ian<fNofMaps; ian++) {
 916     for (Int_t itb=0; itb<fMaxNofSamples; itb++){
 917       Int_t signal=(Int_t)(fHitMap2->GetSignal(ian,itb));
 918       Int_t signalc = Convert10to8(signal);
 919       Int_t signale = Convert8to10(signalc);
 920       AddDigit(ian,itb,signalc,signale);
 921     }
 922   }
 923 }
 924 //______________________________________________________________________
 925 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
 926     // Creates histograms of maps for debugging
 927     Int_t i;
 928
 929     fHis=new TObjArray(fNofMaps);
 930     for (i=0;i<fNofMaps;i++) {
 931         TString sddName("sdd_");
 932         Char_t candNum[4];
 933         sprintf(candNum,"%d",i+1);
 934         sddName.Append(candNum);
 935         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
 936                              0.,(Float_t) scale*fMaxNofSamples), i);
 937     } // end for i
 938 }
 939 //______________________________________________________________________
 940 void AliITSsimulationSDD::FillHistograms(){
 941     // fill 1D histograms from map
 942
 943     if (!fHis) return;
 944
 945     for( Int_t i=0; i<fNofMaps; i++) {
 946         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
 947         Int_t nsamples = hist->GetNbinsX();
 948         for( Int_t j=0; j<nsamples; j++) {
 949             Double_t signal=fHitMap2->GetSignal(i,j);
 950             hist->Fill((Float_t)j,signal);
 951         } // end for j
 952     } // end for i
 953 }
 954 //______________________________________________________________________
 955 void AliITSsimulationSDD::ResetHistograms(){
 956     // Reset histograms for this detector
 957     Int_t i;
 958
 959     for (i=0;i<fNofMaps;i++ ) {
 960         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
 961     } // end for i
 962 }
 963 //______________________________________________________________________
 964 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
 965     // Fills a histogram from a give anode.
 966
 967     if (!fHis) return 0;
 968
 969     if(wing <=0 || wing > 2) {
 970         Warning("GetAnode","Wrong wing number: %d",wing);
 971         return NULL;
 972     } // end if wing <=0 || wing >2
 973     if(anode <=0 || anode > fNofMaps/2) {
 974         Warning("GetAnode","Wrong anode number: %d",anode);
 975         return NULL;
 976     } // end if ampde <=0 || andoe > fNofMaps/2
 977
 978     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
 979     return (TH1F*)(fHis->At(index));
 980 }
 981 //______________________________________________________________________
 982 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
 983     // Writes the histograms to a file
 984
 985     if (!fHis) return;
 986
 987     hfile->cd();
 988     Int_t i;
 989     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
 990     return;
 991 }
 992 //______________________________________________________________________
 993 void AliITSsimulationSDD::WriteSDigits(){
 994     // Fills the Summable digits Tree
 995     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 996
 997     for( Int_t i=0; i<fNofMaps; i++ ) {
 998         if( !fAnodeFire[i] ) continue;
 999         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1000             Double_t sig = fHitMap2->GetSignal( i, j );
1001             if( sig > 0.2 ) {
1002                 Int_t jdx = j*fScaleSize;
1003                 Int_t index = fpList->GetHitIndex( i, j );
1004                 AliITSpListItem pItemTmp2( fModule, index, 0. );
1005                 // put the fScaleSize analog digits in only one
1006                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1007                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1008                     if( pItemTmp == 0 ) continue;
1009                     pItemTmp2.Add( pItemTmp );
1010                 }
1011                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1012                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1013                 aliITS->AddSumDigit( pItemTmp2 );
1014             } // end if (sig > 0.2)
1015         }
1016     }
1017     return;
1018 }
1019 //______________________________________________________________________
1020 void AliITSsimulationSDD::PrintStatus() const {
1021     // Print SDD simulation Parameters
1022
1023     cout << "**************************************************" << endl;
1024     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1025     cout << "**************************************************" << endl;
1026     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1027     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1028     cout << "Number of Anodes used: " << fNofMaps << endl;
1029     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1030     cout << "Scale size factor: " << fScaleSize << endl;
1031     cout << "**************************************************" << endl;
1032 }