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 Int_t power(Int_t b, Int_t e) {
  54     // compute b to the e power, where both b and e are Int_ts.
  55     Int_t power = 1,i;
  56
  57     for(i=0; i<e; i++) power *= b;
  58     return power;
  59 }
  60 //______________________________________________________________________
  61 void FastFourierTransform(AliITSetfSDD *alisddetf,Double_t *real,
  62                           Double_t *imag,Int_t direction) {
  63     // Do a Fast Fourier Transform
  64
  65     Int_t samples = alisddetf->GetSamples();
  66     Int_t l = (Int_t) ((log((Float_t) samples)/log(2.))+0.5);
  67     Int_t m1 = samples;
  68     Int_t m  = samples/2;
  69     Int_t m2 = samples/m1;
  70     Int_t i,j,k;
  71     for(i=1; i<=l; i++) {
  72         for(j=0; j<samples; j += m1) {
  73             Int_t p = 0;
  74             for(k=j; k<= j+m-1; k++) {
  75                 Double_t wsr = alisddetf->GetWeightReal(p);
  76                 Double_t wsi = alisddetf->GetWeightImag(p);
  77                 if(direction == -1) wsi = -wsi;
  78                 Double_t xr = *(real+k+m);
  79                 Double_t xi = *(imag+k+m);
  80                 *(real+k+m) = wsr*(*(real+k)-xr) - wsi*(*(imag+k)-xi);
  81                 *(imag+k+m) = wsr*(*(imag+k)-xi) + wsi*(*(real+k)-xr);
  82                 *(real+k) += xr;
  83                 *(imag+k) += xi;
  84                 p += m2;
  85             } // end for k
  86         } // end for j
  87         m1 = m;
  88         m /= 2;
  89         m2 += m2;
  90     } // end for i
  91
  92     for(j=0; j<samples; j++) {
  93         Int_t j1 = j;
  94         Int_t p = 0;
  95         Int_t i1;
  96         for(i1=1; i1<=l; i1++) {
  97             Int_t j2 = j1;
  98             j1 /= 2;
  99             p = p + p + j2 - j1 - j1;
 100         } // end for i1
 101         if(p >= j) {
 102             Double_t xr = *(real+j);
 103             Double_t xi = *(imag+j);
 104             *(real+j) = *(real+p);
 105             *(imag+j) = *(imag+p);
 106             *(real+p) = xr;
 107             *(imag+p) = xi;
 108         } // end if p>=j
 109     } // end for j
 110     if(direction == -1) {
 111         for(i=0; i<samples; i++) {
 112             *(real+i) /= samples;
 113             *(imag+i) /= samples;
 114         } // end for i
 115     } // end if direction == -1
 116     return;
 117 }
 118 //______________________________________________________________________
 119 AliITSsimulationSDD::AliITSsimulationSDD():
 120 AliITSsimulation(),
 121 fITS(0),
 122 fHitMap2(0),
 123 fHitSigMap2(0),
 124 fHitNoiMap2(0),
 125 fStream(0),
 126 fElectronics(0),
 127 fInZR(0),
 128 fInZI(0),
 129 fOutZR(0),
 130 fOutZI(0),
 131 fAnodeFire(0),
 132 fHis(0),
 133 fD(),
 134 fT1(),
 135 fT2(),
 136 fTol(),
 137 fTreeB(0),
 138 fParam(0),
 139 fFileName(),
 140 fFlag(kFALSE),
 141 fCheckNoise(kFALSE),
 142 fCrosstalkFlag(kFALSE),
 143 fDoFFT(1),
 144 fNofMaps(0),
 145 fMaxNofSamples(0),
 146 fScaleSize(0){
 147     // Default constructor
 148     SetScaleFourier();
 149     SetPerpendTracksFlag();
 150     SetCrosstalkFlag();
 151     SetDoFFT();
 152     SetCheckNoise();
 153 }
 154 //______________________________________________________________________
 155 AliITSsimulationSDD::AliITSsimulationSDD(const AliITSsimulationSDD &source) :
 156     AliITSsimulation(source),
 157 fITS(source.fITS),
 158 fHitMap2(source.fHitMap2),
 159 fHitSigMap2(source.fHitSigMap2),
 160 fHitNoiMap2(source.fHitNoiMap2),
 161 fStream(source.fStream),
 162 fElectronics(source.fElectronics),
 163 fInZR(source.fInZR),
 164 fInZI(source.fInZI),
 165 fOutZR(source.fOutZR),
 166 fOutZI(source.fOutZI),
 167 fAnodeFire(source.fAnodeFire),
 168 fHis(source.fHis),
 169 fD(source.fD),
 170 fT1(source.fT1),
 171 fT2(source.fT2),
 172 fTol(source.fTol),
 173 fTreeB(source.fTreeB),
 174 fParam(source.fParam),
 175 fFileName(source.fFileName),
 176 fFlag(source.fFlag),
 177 fCheckNoise(source.fCheckNoise),
 178 fCrosstalkFlag(source.fCrosstalkFlag),
 179 fDoFFT(source.fDoFFT),
 180 fNofMaps(source.fNofMaps),
 181 fMaxNofSamples(source.fMaxNofSamples),
 182 fScaleSize(source.fScaleSize){
 183     // Copy constructor to satify Coding roules only.
 184
 185 }
 186 //______________________________________________________________________
 187 AliITSsimulationSDD& AliITSsimulationSDD::operator=(const AliITSsimulationSDD &src){
 188     // Assignment operator to satify Coding roules only.
 189
 190     if(this==&src) return *this;
 191     Error("AliITSsimulationSDD","Not allowed to make a = with "
 192           "AliITSsimulationSDD Using default creater instead");
 193     return *this ;
 194 }
 195 //______________________________________________________________________
 196 AliITSsimulation& AliITSsimulationSDD::operator=(const AliITSsimulation &src){
 197     // Assignment operator to satify Coding roules only.
 198
 199     if(this==&src) return *this;
 200     Error("AliITSsimulationSSD","Not allowed to make a = with "
 201           "AliITSsimulationSDD Using default creater instead");
 202     return *this ;
 203 }
 204
 205 //______________________________________________________________________
 206 AliITSsimulationSDD::AliITSsimulationSDD(AliITSDetTypeSim* dettyp):
 207 AliITSsimulation(dettyp),
 208 fITS(0),
 209 fHitMap2(0),
 210 fHitSigMap2(0),
 211 fHitNoiMap2(0),
 212 fStream(0),
 213 fElectronics(0),
 214 fInZR(0),
 215 fInZI(0),
 216 fOutZR(0),
 217 fOutZI(0),
 218 fAnodeFire(0),
 219 fHis(0),
 220 fD(),
 221 fT1(),
 222 fT2(),
 223 fTol(),
 224 fTreeB(0),
 225 fParam(),
 226 fFileName(),
 227 fFlag(kFALSE),
 228 fCheckNoise(kFALSE),
 229 fCrosstalkFlag(kFALSE),
 230 fDoFFT(1),
 231 fNofMaps(0),
 232 fMaxNofSamples(0),
 233 fScaleSize(0){
 234     // Default Constructor
 235   Init();
 236 }
 237 //______________________________________________________________________
 238 void AliITSsimulationSDD::Init(){
 239     // Standard Constructor
 240
 241     SetScaleFourier();
 242     SetPerpendTracksFlag();
 243     SetCrosstalkFlag();
 244     SetDoFFT();
 245     SetCheckNoise();
 246
 247     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 248
 249     AliITSresponseSDD* res = (AliITSresponseSDD*)fDetType->GetResponse(1);
 250     fpList = new AliITSpList( seg->Npz(),
 251                               fScaleSize*seg->Npx() );
 252     fHitSigMap2 = new AliITSMapA2(seg,fScaleSize,1);
 253     fHitNoiMap2 = new AliITSMapA2(seg,fScaleSize,1);
 254     fHitMap2 = fHitSigMap2;
 255
 256     fNofMaps = seg->Npz();
 257     fMaxNofSamples = seg->Npx();
 258     fAnodeFire = new Bool_t [fNofMaps];
 259
 260     Float_t sddLength = seg->Dx();
 261     Float_t sddWidth  = seg->Dz();
 262
 263     Int_t dummy        = 0;
 264     Float_t anodePitch = seg->Dpz(dummy);
 265     Double_t timeStep  = (Double_t)seg->Dpx(dummy);
 266     Float_t driftSpeed = res->DriftSpeed();
 267
 268     if(anodePitch*(fNofMaps/2) > sddWidth) {
 269         Warning("AliITSsimulationSDD",
 270                 "Too many anodes %d or too big pitch %f \n",
 271                 fNofMaps/2,anodePitch);
 272     } // end if
 273
 274     if(timeStep*fMaxNofSamples < sddLength/driftSpeed) {
 275         Error("AliITSsimulationSDD",
 276               "Time Interval > Allowed Time Interval: exit\n");
 277         return;
 278     } // end if
 279
 280     fElectronics = new AliITSetfSDD(timeStep/fScaleSize,
 281                                     res->Electronics());
 282
 283     char opt1[20], opt2[20];
 284     res->ParamOptions(opt1,opt2);
 285     fParam = opt2;
 286
 287     const char *kopt=res->ZeroSuppOption();
 288     fD.Set(fNofMaps);
 289     fT1.Set(fNofMaps);
 290     fT2.Set(fNofMaps);
 291     fTol.Set(fNofMaps);
 292
 293     Bool_t write = res->OutputOption();
 294     if(write && strstr(kopt,"2D")) MakeTreeB();
 295
 296     fITS       = (AliITS*)gAlice->GetModule("ITS");
 297     Int_t size = fNofMaps*fMaxNofSamples;
 298     fStream    = new AliITSInStream(size);
 299
 300     fInZR  = new Double_t [fScaleSize*fMaxNofSamples];
 301     fInZI  = new Double_t [fScaleSize*fMaxNofSamples];
 302     fOutZR = new Double_t [fScaleSize*fMaxNofSamples];
 303     fOutZI = new Double_t [fScaleSize*fMaxNofSamples];
 304 }
 305 //______________________________________________________________________
 306 AliITSsimulationSDD::~AliITSsimulationSDD() {
 307     // destructor
 308
 309     //    delete fpList;
 310     delete fHitSigMap2;
 311     delete fHitNoiMap2;
 312     delete fStream;
 313     delete fElectronics;
 314
 315     fITS = 0;
 316
 317     if (fHis) {
 318         fHis->Delete();
 319         delete fHis;
 320     } // end if fHis
 321     if(fTreeB) delete fTreeB;
 322     if(fInZR)  delete [] fInZR;
 323     if(fInZI)  delete [] fInZI;
 324     if(fOutZR) delete [] fOutZR;
 325     if(fOutZI) delete [] fOutZI;
 326     if(fAnodeFire) delete [] fAnodeFire;
 327 }
 328 //______________________________________________________________________
 329 void AliITSsimulationSDD::InitSimulationModule( Int_t module, Int_t event ) {
 330     // create maps to build the lists of tracks for each summable digit
 331     fModule = module;
 332     fEvent  = event;
 333     ClearMaps();
 334     memset(fAnodeFire,0,sizeof(Bool_t)*fNofMaps);
 335 }
 336 //______________________________________________________________________
 337 void AliITSsimulationSDD::ClearMaps() {
 338     // clear maps
 339     fpList->ClearMap();
 340     fHitSigMap2->ClearMap();
 341     fHitNoiMap2->ClearMap();
 342 }
 343 //______________________________________________________________________
 344 void AliITSsimulationSDD::SDigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 345     // digitize module using the "slow" detector simulator creating
 346     // summable digits.
 347
 348     TObjArray *fHits = mod->GetHits();
 349     Int_t nhits      = fHits->GetEntriesFast();
 350     if( !nhits ) return;
 351
 352     InitSimulationModule( md, ev );
 353     HitsToAnalogDigits( mod );
 354     ChargeToSignal( fModule,kFALSE ); // - Process signal without add noise
 355     fHitMap2 = fHitNoiMap2;   // - Swap to noise map
 356     ChargeToSignal( fModule,kTRUE );  // - Process only noise
 357     fHitMap2 = fHitSigMap2;   // - Return to signal map
 358     WriteSDigits();
 359     ClearMaps();
 360 }
 361 //______________________________________________________________________
 362 Bool_t AliITSsimulationSDD::AddSDigitsToModule(TClonesArray *pItemArray,
 363                                                Int_t mask ) {
 364     // Add Summable digits to module maps.
 365    AliITSresponseSDD* res = (AliITSresponseSDD*)fDetType->GetResponse(1);
 366     Int_t    nItems = pItemArray->GetEntries();
 367     Double_t maxadc = res->MaxAdc();
 368     Bool_t sig = kFALSE;
 369
 370     // cout << "Adding "<< nItems <<" SDigits to module " << fModule << endl;
 371     for( Int_t i=0; i<nItems; i++ ) {
 372         AliITSpListItem * pItem = (AliITSpListItem *)(pItemArray->At( i ));
 373         if( pItem->GetModule() != fModule ) {
 374             Error( "AliITSsimulationSDD","Error reading, SDigits module "
 375                    "%d != current module %d: exit",
 376                    pItem->GetModule(), fModule );
 377             return sig;
 378         } // end if
 379
 380         if(pItem->GetSignal()>0.0 ) sig = kTRUE;
 381
 382         fpList->AddItemTo( mask, pItem ); // Add SignalAfterElect + noise
 383         AliITSpListItem * pItem2 = fpList->GetpListItem( pItem->GetIndex() );
 384         Double_t sigAE = pItem2->GetSignalAfterElect();
 385         if( sigAE >= maxadc ) sigAE = maxadc-1; // avoid overflow signal
 386         Int_t ia;
 387         Int_t it;
 388         fpList->GetMapIndex( pItem->GetIndex(), ia, it );
 389         fHitMap2->SetHit( ia, it, sigAE );
 390         fAnodeFire[ia] = kTRUE;
 391     }
 392     return sig;
 393 }
 394 //______________________________________________________________________
 395 void AliITSsimulationSDD::FinishSDigitiseModule() {
 396     // digitize module using the "slow" detector simulator from
 397     // the sum of summable digits.
 398     FinishDigits() ;
 399     ClearMaps();
 400 }
 401 //______________________________________________________________________
 402 void AliITSsimulationSDD::DigitiseModule(AliITSmodule *mod,Int_t md,Int_t ev){
 403     // create maps to build the lists of tracks for each digit
 404
 405     TObjArray *fHits = mod->GetHits();
 406     Int_t nhits      = fHits->GetEntriesFast();
 407
 408     InitSimulationModule( md, ev );
 409
 410     if( !nhits && fCheckNoise ) {
 411         ChargeToSignal( fModule,kTRUE );  // process noise
 412         GetNoise();
 413         ClearMaps();
 414         return;
 415     } else
 416         if( !nhits ) return;
 417
 418     HitsToAnalogDigits( mod );
 419     ChargeToSignal( fModule,kTRUE );  // process signal + noise
 420
 421     for( Int_t i=0; i<fNofMaps; i++ ) {
 422         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
 423             Int_t jdx = j*fScaleSize;
 424             Int_t index = fpList->GetHitIndex( i, j );
 425             AliITSpListItem pItemTmp2( fModule, index, 0. );
 426             // put the fScaleSize analog digits in only one
 427             for( Int_t ik=0; ik<fScaleSize; ik++ ) {
 428                 AliITSpListItem *pItemTmp = fpList->GetpListItem( i, jdx+ik );
 429                 if( pItemTmp == 0 ) continue;
 430                 pItemTmp2.Add( pItemTmp );
 431             }
 432             fpList->DeleteHit( i, j );
 433             fpList->AddItemTo( 0, &pItemTmp2 );
 434         }
 435     }
 436     FinishDigits();
 437     ClearMaps();
 438 }
 439 //______________________________________________________________________
 440 void AliITSsimulationSDD::FinishDigits() {
 441     // introduce the electronics effects and do zero-suppression if required
 442
 443     ApplyDeadChannels(fModule);
 444     if( fCrosstalkFlag ) ApplyCrosstalk(fModule);
 445
 446     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 447     const char *kopt = res->GetZeroSuppOption();
 448     ZeroSuppression( kopt );
 449 }
 450 //______________________________________________________________________
 451 void AliITSsimulationSDD::HitsToAnalogDigits( AliITSmodule *mod ) {
 452     // create maps to build the lists of tracks for each digit
 453
 454   AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 455   AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 456   TObjArray *hits     = mod->GetHits();
 457     Int_t      nhits    = hits->GetEntriesFast();
 458
 459     //    Int_t      arg[6]   = {0,0,0,0,0,0};
 460     Int_t     dummy      = 0;
 461     Int_t     nofAnodes  = fNofMaps/2;
 462     Double_t  sddLength  = seg->Dx();
 463     Double_t  sddWidth   = seg->Dz();
 464     Double_t  anodePitch = seg->Dpz(dummy);
 465     Double_t  timeStep   = seg->Dpx(dummy);
 466     Double_t  driftSpeed = res->GetDriftSpeed();
 467     //Float_t   maxadc     = res->GetMaxAdc();
 468     //Float_t   topValue   = res->GetDynamicRange();
 469     Double_t  norm       = res->GetMaxAdc()/res->GetDynamicRange(); //   maxadc/topValue;
 470     Double_t  cHloss     = res->GetChargeLoss();
 471     Float_t   dfCoeff, s1; res->DiffCoeff(dfCoeff,s1); // Signal 2d Shape
 472     Double_t  eVpairs    = res->GetGeVToCharge()*1.0E9; // 3.6 eV by def.
 473     Double_t  nsigma     = res->GetNSigmaIntegration(); //
 474     Int_t     nlookups   = res->GetGausNLookUp();       //
 475     Float_t   jitter     = res->GetJitterError(); //
 476
 477     // Piergiorgio's part (apart for few variables which I made float
 478     // when i thought that can be done
 479     // Fill detector maps with GEANT hits
 480     // loop over hits in the module
 481
 482     const Float_t kconv = 1.0e+6;  // GeV->KeV
 483     Int_t     itrack      = 0;
 484     Int_t     hitDetector; // detector number (lay,lad,hitDetector)
 485     Int_t     iWing;       // which detector wing/side.
 486     Int_t     detector;    // 2*(detector-1)+iWing
 487     Int_t     ii,kk,ka,kt; // loop indexs
 488     Int_t     ia,it,index; // sub-pixel integration indexies
 489     Int_t     iAnode;      // anode number.
 490     Int_t     timeSample;  // time buckett.
 491     Int_t     anodeWindow; // anode direction charge integration width
 492     Int_t     timeWindow;  // time direction charge integration width
 493     Int_t     jamin,jamax; // anode charge integration window
 494     Int_t     jtmin,jtmax; // time charge integration window
 495     Int_t     ndiv;        // Anode window division factor.
 496     Int_t     nsplit;      // the number of splits in anode and time windows==1.
 497     Int_t     nOfSplits;   // number of times track length is split into
 498     Float_t   nOfSplitsF;  // Floating point version of nOfSplits.
 499     Float_t   kkF;         // Floating point version of loop index kk.
 500     Double_t  pathInSDD; // Track length in SDD.
 501     Double_t  drPath; // average position of track in detector. in microns
 502     Double_t  drTime; // Drift time
 503     Double_t  nmul;   // drift time window multiplication factor.
 504     Double_t  avDrft;  // x position of path length segment in cm.
 505     Double_t  avAnode; // Anode for path length segment in Anode number (float)
 506     Double_t  xAnode;  // Floating point anode number.
 507     Double_t  driftPath; // avDrft in microns.
 508     Double_t  width;     // width of signal at anodes.
 509     Double_t  depEnergy; // Energy deposited in this GEANT step.
 510     Double_t  xL[3],dxL[3]; // local hit coordinates and diff.
 511     Double_t  sigA; // sigma of signal at anode.
 512     Double_t  sigT; // sigma in time/drift direction for track segment
 513     Double_t  aStep,aConst; // sub-pixel size and offset anode
 514     Double_t  tStep,tConst; // sub-pixel size and offset time
 515     Double_t  amplitude; // signal amplitude for track segment in nanoAmpere
 516     Double_t  chargeloss; // charge loss for track segment.
 517     Double_t  anodeAmplitude; // signal amplitude in anode direction
 518     Double_t  aExpo;          // exponent of Gaussian anode direction
 519     Double_t  timeAmplitude;  // signal amplitude in time direction
 520     Double_t  tExpo;          // exponent of Gaussian time direction
 521     //  Double_t tof;            // Time of flight in ns of this step.
 522
 523     for(ii=0; ii<nhits; ii++) {
 524         if(!mod->LineSegmentL(ii,xL[0],dxL[0],xL[1],dxL[1],xL[2],dxL[2],
 525                               depEnergy,itrack)) continue;
 526         xL[0] += 0.0001*gRandom->Gaus( 0, jitter ); //
 527         depEnergy  *= kconv;
 528         hitDetector = mod->GetDet();
 529         //tof         = 1.E+09*(mod->GetHit(ii)->GetTOF()); // tof in ns.
 530         //if(tof>sddLength/driftSpeed) continue; // hit happed too late.
 531
 532         // scale path to simulate a perpendicular track
 533         // continue if the particle did not lose energy
 534         // passing through detector
 535         if (!depEnergy) {
 536             AliDebug(1,
 537              Form("fTrack = %d hit=%d module=%d This particle has passed without losing energy!",
 538                   itrack,ii,mod->GetIndex()));
 539             continue;
 540         } // end if !depEnergy
 541
 542         pathInSDD = TMath::Sqrt(dxL[0]*dxL[0]+dxL[1]*dxL[1]+dxL[2]*dxL[2]);
 543
 544         if (fFlag && pathInSDD) { depEnergy *= (0.03/pathInSDD); }
 545         drPath = 10000.*(dxL[0]+2.*xL[0])*0.5;
 546         if(drPath < 0) drPath = -drPath;
 547         drPath = sddLength-drPath;
 548         if(drPath < 0) {
 549             AliDebug(1, // this should be fixed at geometry level
 550               Form("negative drift path drPath=%e sddLength=%e dxL[0]=%e xL[0]=%e",
 551                    drPath,sddLength,dxL[0],xL[0]));
 552             continue;
 553         } // end if drPath < 0
 554
 555         // Compute number of segments to brake step path into
 556         drTime = drPath/driftSpeed;  //   Drift Time
 557         sigA   = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);// Sigma along the anodes
 558         // calcuate the number of time the path length should be split into.
 559         nOfSplits = (Int_t) (1. + 10000.*pathInSDD/sigA);
 560         if(fFlag) nOfSplits = 1;
 561
 562         // loop over path segments, init. some variables.
 563         depEnergy /= nOfSplits;
 564         nOfSplitsF = (Float_t) nOfSplits;
 565         for(kk=0;kk<nOfSplits;kk++) { // loop over path segments
 566             kkF       = (Float_t) kk + 0.5;
 567             avDrft    = xL[0]+dxL[0]*kkF/nOfSplitsF;
 568             avAnode   = xL[2]+dxL[2]*kkF/nOfSplitsF;
 569             driftPath = 10000.*avDrft;
 570
 571             iWing = 2;  // Assume wing is 2
 572             if(driftPath < 0) { // if wing is not 2 it is 1.
 573                 iWing     = 1;
 574                 driftPath = -driftPath;
 575             } // end if driftPath < 0
 576             driftPath = sddLength-driftPath;
 577             detector  = 2*(hitDetector-1) + iWing;
 578             if(driftPath < 0) {
 579                 AliDebug(1, // this should be fixed at geometry level
 580                  Form("negative drift path driftPath=%e sddLength=%e avDrft=%e dxL[0]=%e xL[0]=%e",
 581                       driftPath,sddLength,avDrft,dxL[0],xL[0]));
 582                 continue;
 583             } // end if driftPath < 0
 584
 585             //   Drift Time
 586             drTime     = driftPath/driftSpeed; // drift time for segment.
 587             timeSample = (Int_t) (fScaleSize*drTime/timeStep + 1);
 588             // compute time Sample including tof information. The tof only
 589             // effects the time of the signal is recoreded and not the
 590             // the defusion.
 591             // timeSample = (Int_t) (fScaleSize*(drTime+tof)/timeStep + 1);
 592             if(timeSample > fScaleSize*fMaxNofSamples) {
 593                 Warning("HitsToAnalogDigits","Wrong Time Sample: %e",
 594                         timeSample);
 595                 continue;
 596             } // end if timeSample > fScaleSize*fMaxNoofSamples
 597
 598             //   Anode
 599             xAnode = 10000.*(avAnode)/anodePitch + nofAnodes/2;  // +1?
 600             if(xAnode*anodePitch > sddWidth || xAnode*anodePitch < 0.)
 601                 Warning("HitsToAnalogDigits",
 602                         "Exceedubg sddWidth=%e Z = %e",
 603                         sddWidth,xAnode*anodePitch);
 604             iAnode = (Int_t) (1.+xAnode); // xAnode?
 605             if(iAnode < 1 || iAnode > nofAnodes) {
 606                 Warning("HitToAnalogDigits","Wrong iAnode: 1<%d>%d  (xanode=%e)",
 607                         iAnode,nofAnodes, xAnode);
 608                 continue;
 609             } // end if iAnode < 1 || iAnode > nofAnodes
 610
 611             // store straight away the particle position in the array
 612             // of particles and take idhit=ii only when part is entering (this
 613             // requires FillModules() in the macro for analysis) :
 614
 615             // Sigma along the anodes for track segment.
 616             sigA       = TMath::Sqrt(2.*dfCoeff*drTime+s1*s1);
 617             sigT       = sigA/driftSpeed;
 618             // Peak amplitude in nanoAmpere
 619             amplitude  = fScaleSize*160.*depEnergy/
 620                 (timeStep*eVpairs*2.*acos(-1.)*sigT*sigA);
 621             amplitude *= timeStep/25.; // WARNING!!!!! Amplitude scaling to
 622             // account for clock variations
 623             // (reference value: 40 MHz)
 624             chargeloss = 1.-cHloss*driftPath/1000;
 625             amplitude *= chargeloss;
 626             width  = 2.*nsigma/(nlookups-1);
 627             // Spread the charge
 628             // Pixel index
 629             ndiv = 2;
 630             nmul = 3.;
 631             if(drTime > 1200.) {
 632                 ndiv = 4;
 633                 nmul = 1.5;
 634             } // end if drTime > 1200.
 635             // Sub-pixel index
 636             nsplit = 4; // hard-wired //nsplit=4;nsplit = (nsplit+1)/2*2;
 637             // Sub-pixel size see computation of aExpo and tExpo.
 638             aStep  = anodePitch/(nsplit*fScaleSize*sigA);
 639             aConst = xAnode*anodePitch/sigA;
 640             tStep  = timeStep/(nsplit*fScaleSize*sigT);
 641             tConst = drTime/sigT;
 642             // Define SDD window corresponding to the hit
 643             anodeWindow = (Int_t)(fScaleSize*nsigma*sigA/anodePitch+1);
 644             timeWindow  = (Int_t) (fScaleSize*nsigma*sigT/timeStep+1.);
 645             jamin = (iAnode - anodeWindow/ndiv - 1)*fScaleSize*nsplit +1;
 646             jamax = (iAnode + anodeWindow/ndiv)*fScaleSize*nsplit;
 647             if(jamin <= 0) jamin = 1;
 648             if(jamax > fScaleSize*nofAnodes*nsplit)
 649                 jamax = fScaleSize*nofAnodes*nsplit;
 650             // jtmin and jtmax are Hard-wired
 651             jtmin = (Int_t)(timeSample-timeWindow*nmul-1)*nsplit+1;
 652             jtmax = (Int_t)(timeSample+timeWindow*nmul)*nsplit;
 653             if(jtmin <= 0) jtmin = 1;
 654             if(jtmax > fScaleSize*fMaxNofSamples*nsplit)
 655                 jtmax = fScaleSize*fMaxNofSamples*nsplit;
 656             // Spread the charge in the anode-time window
 657             for(ka=jamin; ka <=jamax; ka++) {
 658                 ia = (ka-1)/(fScaleSize*nsplit) + 1;
 659                 if(ia <= 0) {
 660                     Warning("HitsToAnalogDigits","ia < 1: ");
 661                     continue;
 662                 } // end if
 663                 if(ia > nofAnodes) ia = nofAnodes;
 664                 aExpo     = (aStep*(ka-0.5)-aConst);
 665                 if(TMath::Abs(aExpo) > nsigma)  anodeAmplitude = 0.;
 666                 else {
 667                     dummy          = (Int_t) ((aExpo+nsigma)/width);
 668                     anodeAmplitude = amplitude*res->GetGausLookUp(dummy);
 669                 } // end if TMath::Abs(aEspo) > nsigma
 670                 // index starts from 0
 671                 index = ((detector+1)%2)*nofAnodes+ia-1;
 672                 if(anodeAmplitude) for(kt=jtmin; kt<=jtmax; kt++) {
 673                     it = (kt-1)/nsplit+1;  // it starts from 1
 674                     if(it<=0){
 675                         Warning("HitsToAnalogDigits","it < 1:");
 676                         continue;
 677                     } // end if
 678                     if(it>fScaleSize*fMaxNofSamples)
 679                         it = fScaleSize*fMaxNofSamples;
 680                     tExpo    = (tStep*(kt-0.5)-tConst);
 681                     if(TMath::Abs(tExpo) > nsigma) timeAmplitude = 0.;
 682                     else {
 683                         dummy         = (Int_t) ((tExpo+nsigma)/width);
 684                         timeAmplitude = anodeAmplitude*
 685                             res->GetGausLookUp(dummy);
 686                     } // end if TMath::Abs(tExpo) > nsigma
 687                     // build the list of Sdigits for this module
 688                     //                    arg[0]     = index;
 689                     //                    arg[1]     = it;
 690                     //                    arg[2]     = itrack; // track number
 691                     //                    arg[3]     = ii-1; // hit number.
 692                     timeAmplitude *= norm;
 693                     timeAmplitude *= 10;
 694                     //         ListOfFiredCells(arg,timeAmplitude,alst,padr);
 695                     Double_t charge = timeAmplitude;
 696                     charge += fHitMap2->GetSignal(index,it-1);
 697                     fHitMap2->SetHit(index, it-1, charge);
 698                     fpList->AddSignal(index,it-1,itrack,ii-1,
 699                                       mod->GetIndex(),timeAmplitude);
 700                     fAnodeFire[index] = kTRUE;
 701                 } // end if anodeAmplitude and loop over time in window
 702             } // loop over anodes in window
 703         } // end loop over "sub-hits"
 704     } // end loop over hits
 705 }
 706 /*
 707 //______________________________________________________________________
 708 void AliITSsimulationSDD::ListOfFiredCells(Int_t *arg,Double_t timeAmplitude,
 709                                            TObjArray *alist,TClonesArray *padr){
 710     // Returns the list of "fired" cells.
 711
 712     Int_t index     = arg[0];
 713     Int_t ik        = arg[1];
 714     Int_t idtrack   = arg[2];
 715     Int_t idhit     = arg[3];
 716     Int_t counter   = arg[4];
 717     Int_t countadr  = arg[5];
 718     Double_t charge = timeAmplitude;
 719     charge += fHitMap2->GetSignal(index,ik-1);
 720     fHitMap2->SetHit(index, ik-1, charge);
 721
 722     Int_t digits[3];
 723     Int_t it = (Int_t)((ik-1)/fScaleSize);
 724     digits[0] = index;
 725     digits[1] = it;
 726     digits[2] = (Int_t)timeAmplitude;
 727     Float_t phys;
 728     if (idtrack >= 0) phys = (Float_t)timeAmplitude;
 729     else phys = 0;
 730
 731     Double_t cellcharge = 0.;
 732     AliITSTransientDigit* pdigit;
 733     // build the list of fired cells and update the info
 734     if (!fHitMap1->TestHit(index, it)) {
 735         new((*padr)[countadr++]) TVector(3);
 736         TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
 737         trinfo(0) = (Float_t)idtrack;
 738         trinfo(1) = (Float_t)idhit;
 739         trinfo(2) = (Float_t)timeAmplitude;
 740
 741         alist->AddAtAndExpand(new AliITSTransientDigit(phys,digits),counter);
 742         fHitMap1->SetHit(index, it, counter);
 743         counter++;
 744         pdigit=(AliITSTransientDigit*)alist->At(alist->GetLast());
 745         // list of tracks
 746         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
 747         trlist->Add(&trinfo);
 748     } else {
 749         pdigit = (AliITSTransientDigit*) fHitMap1->GetHit(index, it);
 750         for(Int_t kk=0;kk<fScaleSize;kk++) {
 751             cellcharge += fHitMap2->GetSignal(index,fScaleSize*it+kk);
 752         }  // end for kk
 753         // update charge
 754         (*pdigit).fSignal = (Int_t)cellcharge;
 755         (*pdigit).fPhysics += phys;
 756         // update list of tracks
 757         TObjArray* trlist = (TObjArray*)pdigit->TrackList();
 758         Int_t lastentry = trlist->GetLast();
 759         TVector *ptrkp = (TVector*)trlist->At(lastentry);
 760         TVector &trinfo = *ptrkp;
 761         Int_t lasttrack = Int_t(trinfo(0));
 762         Float_t lastcharge=(trinfo(2));
 763         if (lasttrack==idtrack ) {
 764             lastcharge += (Float_t)timeAmplitude;
 765             trlist->RemoveAt(lastentry);
 766             trinfo(0) = lasttrack;
 767             trinfo(1) = idhit;
 768             trinfo(2) = lastcharge;
 769             trlist->AddAt(&trinfo,lastentry);
 770         } else {
 771             new((*padr)[countadr++]) TVector(3);
 772             TVector &trinfo=*((TVector*) (*padr)[countadr-1]);
 773             trinfo(0) = (Float_t)idtrack;
 774             trinfo(1) = (Float_t)idhit;
 775             trinfo(2) = (Float_t)timeAmplitude;
 776             trlist->Add(&trinfo);
 777         } // end if lasttrack==idtrack
 778
 779         if(AliDebugLevel()){
 780             // check the track list - debugging
 781             Int_t trk[20], htrk[20];
 782             Float_t chtrk[20];
 783             Int_t nptracks = trlist->GetEntriesFast();
 784             if (nptracks > 2) {
 785                 Int_t tr;
 786                 for (tr=0;tr<nptracks;tr++) {
 787                     TVector *pptrkp = (TVector*)trlist->At(tr);
 788                     TVector &pptrk  = *pptrkp;
 789                     trk[tr]   = Int_t(pptrk(0));
 790                     htrk[tr]  = Int_t(pptrk(1));
 791                     chtrk[tr] = (pptrk(2));
 792                     cout << "nptracks "<<nptracks << endl;
 793                 } // end for tr
 794             } // end if nptracks
 795         } // end if AliDebugLevel()
 796     } //  end if pdigit
 797
 798     // update counter and countadr for next call.
 799     arg[4] = counter;
 800     arg[5] = countadr;
 801 }
 802 */
 803 //____________________________________________
 804 void AliITSsimulationSDD::AddDigit( Int_t i, Int_t j, Int_t signal ) {
 805     // Adds a Digit.
 806     Int_t size = AliITSdigit::GetNTracks();
 807
 808     Int_t digits[3];
 809     Int_t * tracks = new Int_t[size];
 810     Int_t * hits = new Int_t[size];
 811     Float_t phys;
 812     Float_t * charges = new Float_t[size];
 813
 814     digits[0] = i;
 815     digits[1] = j;
 816     digits[2] = signal;
 817
 818     AliITSpListItem *pItem = fpList->GetpListItem( i, j );
 819     if( pItem == 0 ) {
 820         phys = 0.0;
 821         for( Int_t l=0; l<size; l++ ) {
 822             tracks[l]  = 0;
 823             hits[l]    = 0;
 824             charges[l] = 0.0;
 825         }
 826     } else {
 827         Int_t idtrack =  pItem->GetTrack( 0 );
 828         if( idtrack >= 0 ) phys = pItem->GetSignal();
 829         else phys = 0.0;
 830
 831         for( Int_t l=0; l<size; l++ ) if(l<pItem->GetMaxKept()) {
 832             tracks[l]  = pItem->GetTrack( l );
 833             hits[l]    = pItem->GetHit( l );
 834             charges[l] = pItem->GetSignal( l );
 835         }else{
 836             tracks[l]  = -3;
 837             hits[l]    = -1;
 838             charges[l] = 0.0;
 839         }// end for if
 840     }
 841
 842     fITS->AddSimDigit( 1, phys, digits, tracks, hits, charges );
 843     delete [] tracks;
 844     delete [] hits;
 845     delete [] charges;
 846 }
 847 //______________________________________________________________________
 848 void AliITSsimulationSDD::ChargeToSignal(Int_t mod,Bool_t bAddNoise) {
 849     // add baseline, noise, electronics and ADC saturation effects
 850
 851     char opt1[20], opt2[20];
 852     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 853     res->GetParamOptions(opt1,opt2);
 854     Double_t baseline=0;
 855     Double_t noise=0;
 856
 857     Float_t contrib=0;
 858     Int_t i,k,kk;
 859     Float_t maxadc = res->GetMaxAdc();
 860     if(!fDoFFT) {
 861         for (i=0;i<fNofMaps;i++) {
 862             if( !fAnodeFire[i] ) continue;
 863             baseline = res->GetBaseline(i);
 864             noise = res->GetNoise(i);
 865
 866             for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 867                 fInZR[k]  = fHitMap2->GetSignal(i,k);
 868                 if( bAddNoise ) {
 869                     contrib   = (baseline + noise*gRandom->Gaus());
 870                     fInZR[k] += contrib;
 871                 }
 872             } // end for k
 873             for(k=0; k<fMaxNofSamples; k++) {
 874                 Double_t newcont = 0.;
 875                 Double_t maxcont = 0.;
 876                 for(kk=0;kk<fScaleSize;kk++) {
 877                     newcont = fInZR[fScaleSize*k+kk];
 878                     if(newcont > maxcont) maxcont = newcont;
 879                 } // end for kk
 880                 newcont = maxcont;
 881                 if (newcont >= maxadc) newcont = maxadc -1;
 882                 if(newcont >= baseline){
 883                     Warning("","newcont=%d>=baseline=%d",newcont,baseline);
 884                 } // end if
 885                 // back to analog: ?
 886                 fHitMap2->SetHit(i,k,newcont);
 887             }  // end for k
 888         } // end for i loop over anodes
 889         return;
 890     } // end if DoFFT
 891
 892     for (i=0;i<fNofMaps;i++) {
 893         if( !fAnodeFire[i] ) continue;
 894         baseline = res->GetBaseline(i);
 895         noise = res->GetNoise(i);
 896         for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 897             fInZR[k]  = fHitMap2->GetSignal(i,k);
 898             if( bAddNoise ) {
 899                 contrib   = (baseline + noise*gRandom->Gaus());
 900                 fInZR[k] += contrib;
 901             }
 902             fInZI[k]  = 0.;
 903         } // end for k
 904         FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1);
 905         for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 906             Double_t rw = fElectronics->GetTraFunReal(k);
 907             Double_t iw = fElectronics->GetTraFunImag(k);
 908             fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 909             fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 910         } // end for k
 911         FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1);
 912         for(k=0; k<fMaxNofSamples; k++) {
 913             Double_t newcont1 = 0.;
 914             Double_t maxcont1 = 0.;
 915             for(kk=0;kk<fScaleSize;kk++) {
 916                 newcont1 = fOutZR[fScaleSize*k+kk];
 917                 if(newcont1 > maxcont1) maxcont1 = newcont1;
 918             } // end for kk
 919             newcont1 = maxcont1;
 920             if (newcont1 >= maxadc) newcont1 = maxadc -1;
 921             fHitMap2->SetHit(i,k,newcont1);
 922         } // end for k
 923     } // end for i loop over anodes
 924     return;
 925 }
 926 //____________________________________________________________________
 927 void AliITSsimulationSDD::ApplyDeadChannels(Int_t mod) {
 928     // Set dead channel signal to zero
 929     AliITSCalibrationSDD * calibr = (AliITSCalibrationSDD *)GetCalibrationModel(mod);
 930     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 931     // nothing to do
 932     if( calibr->IsDead() ||
 933         ( calibr->GetDeadChips() == 0 &&
 934           calibr->GetDeadChannels() == 0 ) ) return;
 935
 936     // static AliITS *iTS = (AliITS*)gAlice->GetModule( "ITS" );
 937
 938     Int_t fMaxNofSamples = seg->Npx();
 939     // AliITSgeom *geom = iTS->GetITSgeom();
 940     // Int_t firstSDDMod = geom->GetStartDet( 1 );
 941     // loop over wings
 942     for( Int_t j=0; j<2; j++ ) {
 943       // Int_t mod = (fModule-firstSDDMod)*2 + j;
 944       for( Int_t u=0; u<calibr->Chips(); u++ )
 945         for( Int_t v=0; v<calibr->Channels(); v++ ) {
 946           Float_t gain = calibr->Gain(j, u, v );
 947           for( Int_t k=0; k<fMaxNofSamples; k++ ) {
 948             Int_t i = j*calibr->Chips()*calibr->Channels() +
 949               u*calibr->Channels() +
 950               v;
 951             Double_t signal =  gain * fHitMap2->GetSignal( i, k );
 952             fHitMap2->SetHit( i, k, signal );  ///
 953           }
 954         }
 955     }
 956 }
 957 //______________________________________________________________________
 958 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 959     // function add the crosstalk effect to signal
 960     // temporal function, should be checked...!!!
 961     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 962
 963     Int_t fNofMaps = seg->Npz();
 964     Int_t fMaxNofSamples = seg->Npx();
 965
 966     // create and inizialice crosstalk map
 967     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 968     if( ctk == NULL ) {
 969         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
 970         return;
 971     }
 972     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 973     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 974     for( Int_t z=0; z<fNofMaps; z++ ) {
 975       Double_t baseline = calibr->GetBaseline(z);
 976         Bool_t on = kFALSE;
 977         Int_t tstart = 0;
 978         Int_t tstop = 0;
 979         Int_t nTsteps = 0;
 980
 981         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 982             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 983             if( fadc > baseline ) {
 984                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 985                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 986                     if( fadc1 < fadc ) continue;
 987                     on = kTRUE;
 988                     nTsteps = 0;
 989                     tstart = l;
 990                 }
 991                 nTsteps++;
 992             }
 993             else { // end fadc > baseline
 994                 if( on == kTRUE ) {
 995                     if( nTsteps > 2 ) {
 996                         tstop = l;
 997                         // make smooth derivative
 998                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 999                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
1000                         if( ctk == NULL ) {
1001                             Error( "ApplyCrosstalk",
1002                                    "no memory for temporal array: exit \n" );
1003                             return;
1004                         }
1005                         for( Int_t i=tstart; i<tstop; i++ ) {
1006                             if( i > 2 && i < fMaxNofSamples-2 )
1007                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
1008                                     -0.1*fHitMap2->GetSignal( z,i-1 )
1009                                     +0.1*fHitMap2->GetSignal( z,i+1 )
1010                                     +0.2*fHitMap2->GetSignal( z,i+2 );
1011                         }
1012
1013                         // add crosstalk contribution to neibourg anodes
1014                         for( Int_t i=tstart; i<tstop; i++ ) {
1015                             Int_t anode = z - 1;
1016                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
1017                             Float_t ctktmp =  -dev[i1] * 0.25;
1018                             if( anode > 0 ) {
1019                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
1020                             }
1021                             anode = z + 1;
1022                             if( anode < fNofMaps ) {
1023                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
1024                             }
1025                         }
1026                         delete [] dev;
1027
1028                     } // if( nTsteps > 2 )
1029                     on = kFALSE;
1030                 }  // if( on == kTRUE )
1031             }  // else
1032         }
1033     }
1034
1035     for( Int_t a=0; a<fNofMaps; a++ )
1036         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
1037             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
1038             fHitMap2->SetHit( a, t, signal );
1039         }
1040
1041     delete [] ctk;
1042 }
1043
1044 //______________________________________________________________________
1045 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
1046                                            Int_t &th) const{
1047     // Returns the compression alogirthm parameters
1048   db=fD[i];
1049   tl=fT1[i];
1050   th=fT2[i];
1051 }
1052 //______________________________________________________________________
1053 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl) const{
1054     // returns the compression alogirthm parameters
1055
1056     db=fD[i];
1057     tl=fT1[i];
1058
1059 }
1060 //______________________________________________________________________
1061 void AliITSsimulationSDD::SetCompressParam(){
1062     // Sets the compression alogirthm parameters
1063     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1064     for(Int_t ian = 0; ian<fNofMaps;ian++){
1065       fD[ian] = (Int_t)(calibr->GetBaseline(ian));
1066       fT1[ian] = (Int_t)(2.*calibr->GetNoiseAfterElectronics(ian)+0.5);
1067       fT2[ian] = 0;   // used by 2D clustering - not defined yet
1068       fTol[ian] = 0; // used by 2D clustering - not defined yet
1069     }
1070 }
1071 //______________________________________________________________________
1072 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
1073     // To the 10 to 8 bit lossive compression.
1074     // code from Davide C. and Albert W.
1075
1076     if (signal < 128)  return signal;
1077     if (signal < 256)  return (128+((signal-128)>>1));
1078     if (signal < 512)  return (192+((signal-256)>>3));
1079     if (signal < 1024) return (224+((signal-512)>>4));
1080     return 0;
1081 }
1082 /*
1083 //______________________________________________________________________
1084 AliITSMap*   AliITSsimulationSDD::HitMap(Int_t i){
1085     //Return the correct map.
1086
1087     return ((i==0)? fHitMap1 : fHitMap2);
1088 }
1089 */
1090 //______________________________________________________________________
1091 void AliITSsimulationSDD::ZeroSuppression(const char *option) {
1092     // perform the zero suppresion
1093     if (strstr(option,"2D")) {
1094         //Init2D();              // activate if param change module by module
1095         Compress2D();
1096     } else if (strstr(option,"1D")) {
1097         //Init1D();              // activate if param change module by module
1098         Compress1D();
1099     } else StoreAllDigits();
1100 }
1101 //______________________________________________________________________
1102 void AliITSsimulationSDD::Init2D(){
1103     // read in and prepare arrays: fD, fT1, fT2
1104     //                         savemu[nanodes], savesigma[nanodes]
1105     // read baseline and noise from file - either a .root file and in this
1106     // case data should be organised in a tree with one entry for each
1107     // module => reading should be done accordingly
1108     // or a classic file and do smth. like this ( code from Davide C. and
1109     // Albert W.) :
1110     // Read 2D zero-suppression parameters for SDD
1111
1112     if (!strstr(fParam.Data(),"file")) return;
1113
1114     Int_t na,pos,tempTh;
1115     Float_t mu,sigma;
1116     Float_t *savemu    = new Float_t [fNofMaps];
1117     Float_t *savesigma = new Float_t [fNofMaps];
1118     char input[100],basel[100],par[100];
1119     char *filtmp;
1120     Double_t tmp1,tmp2;
1121     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1122
1123     res->Thresholds(tmp1,tmp2);
1124     Int_t minval = static_cast<Int_t>(tmp1);
1125
1126     res->Filenames(input,basel,par);
1127     fFileName = par;
1128     //
1129     filtmp = gSystem->ExpandPathName(fFileName.Data());
1130     FILE *param = fopen(filtmp,"r");
1131     na = 0;
1132
1133     if(param) {
1134         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1135             if (pos != na+1) {
1136                 Error("Init2D","Anode number not in increasing order!",filtmp);
1137                 exit(1);
1138             } // end if pos != na+1
1139             savemu[na] = mu;
1140             savesigma[na] = sigma;
1141             if ((2.*sigma) < mu) {
1142                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1143                 mu = 2.0 * sigma;
1144             } else fD[na] = 0;
1145             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1146             if (tempTh < 0) tempTh=0;
1147             fT1[na] = tempTh;
1148             tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval;
1149             if (tempTh < 0) tempTh=0;
1150             fT2[na] = tempTh;
1151             na++;
1152         } // end while
1153     } else {
1154         Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp);
1155         exit(1);
1156     } // end if(param)
1157
1158     fclose(param);
1159     delete [] filtmp;
1160     delete [] savemu;
1161     delete [] savesigma;
1162 }
1163 //______________________________________________________________________
1164 void AliITSsimulationSDD::Compress2D(){
1165     // simple ITS cluster finder -- online zero-suppression conditions
1166
1167     Int_t db,tl,th;
1168     Double_t tmp1,tmp2;
1169     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1170
1171     res->Thresholds(tmp1,tmp2);
1172     Int_t minval   = static_cast<Int_t>(tmp1);
1173     Bool_t write   = res->OutputOption();
1174     Bool_t do10to8 = res->Do10to8();
1175     Int_t nz, nl, nh, low, i, j;
1176     SetCompressParam();
1177     for (i=0; i<fNofMaps; i++) {
1178         CompressionParam(i,db,tl,th);
1179         nz  = 0;
1180         nl  = 0;
1181         nh  = 0;
1182         low = 0;
1183         for (j=0; j<fMaxNofSamples; j++) {
1184             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1185             signal -= db; // if baseline eq. is done here
1186             if (signal <= 0) {nz++; continue;}
1187             if ((signal - tl) < minval) low++;
1188             if ((signal - th) >= minval) {
1189                 nh++;
1190                 Bool_t cond=kTRUE;
1191                 FindCluster(i,j,signal,minval,cond);
1192                 if(cond && j &&
1193                    ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){
1194                     if(do10to8) signal = Convert10to8(signal);
1195                     AddDigit(i,j,signal);
1196                 } // end if cond&&j&&()
1197             } else if ((signal - tl) >= minval) nl++;
1198         } // end for j loop time samples
1199         if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1200     } //end for i loop anodes
1201
1202     char hname[30];
1203     if (write) {
1204         sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1205         TreeB()->Write(hname);
1206         // reset tree
1207         TreeB()->Reset();
1208     } // end if write
1209 }
1210 //______________________________________________________________________
1211 void  AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
1212                                        Int_t minval,Bool_t &cond){
1213     // Find clusters according to the online 2D zero-suppression algorithm
1214     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1215     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
1216
1217     Bool_t do10to8 = res->Do10to8();
1218     Bool_t high    = kFALSE;
1219
1220     fHitMap2->FlagHit(i,j);
1221     //
1222     //  check the online zero-suppression conditions
1223     //
1224     const Int_t kMaxNeighbours = 4;
1225     Int_t nn;
1226     Int_t dbx,tlx,thx;
1227     Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
1228     seg->Neighbours(i,j,&nn,xList,yList);
1229     Int_t in,ix,iy,qns;
1230     for (in=0; in<nn; in++) {
1231         ix=xList[in];
1232         iy=yList[in];
1233         if (fHitMap2->TestHit(ix,iy)==kUnused) {
1234             CompressionParam(ix,dbx,tlx,thx);
1235             Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy));
1236             qn -= dbx; // if baseline eq. is done here
1237             if ((qn-tlx) < minval) {
1238                 fHitMap2->FlagHit(ix,iy);
1239                 continue;
1240             } else {
1241                 if ((qn - thx) >= minval) high=kTRUE;
1242                 if (cond) {
1243                     if(do10to8) signal = Convert10to8(signal);
1244                     AddDigit(i,j,signal);
1245                 } // end if cond
1246                 if(do10to8) qns = Convert10to8(qn);
1247                 else qns=qn;
1248                 if (!high) AddDigit(ix,iy,qns);
1249                 cond=kFALSE;
1250                 if(!high) fHitMap2->FlagHit(ix,iy);
1251             } // end if qn-tlx < minval
1252         } // end if  TestHit
1253     } // end for in loop over neighbours
1254 }
1255 //______________________________________________________________________
1256 void AliITSsimulationSDD::Init1D(){
1257     // this is just a copy-paste of input taken from 2D algo
1258     // Torino people should give input
1259     // Read 1D zero-suppression parameters for SDD
1260
1261     if (!strstr(fParam.Data(),"file")) return;
1262
1263     Int_t na,pos,tempTh;
1264     Float_t mu,sigma;
1265     Float_t *savemu    = new Float_t [fNofMaps];
1266     Float_t *savesigma = new Float_t [fNofMaps];
1267     char input[100],basel[100],par[100];
1268     char *filtmp;
1269     Double_t tmp1,tmp2;
1270     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1271
1272     res->Thresholds(tmp1,tmp2);
1273     Int_t minval = static_cast<Int_t>(tmp1);
1274
1275     res->Filenames(input,basel,par);
1276     fFileName=par;
1277
1278     //  set first the disable and tol param
1279     SetCompressParam();
1280     //
1281     filtmp = gSystem->ExpandPathName(fFileName.Data());
1282     FILE *param = fopen(filtmp,"r");
1283     na = 0;
1284
1285     if (param) {
1286         fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]);
1287         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1288             if (pos != na+1) {
1289                 Error("Init1D","Anode number not in increasing order!",filtmp);
1290                 exit(1);
1291             } // end if pos != na+1
1292             savemu[na]=mu;
1293             savesigma[na]=sigma;
1294             if ((2.*sigma) < mu) {
1295                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1296                 mu = 2.0 * sigma;
1297             } else fD[na] = 0;
1298             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1299             if (tempTh < 0) tempTh=0;
1300             fT1[na] = tempTh;
1301             na++;
1302         } // end while
1303     } else {
1304         Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp);
1305         exit(1);
1306     } // end if(param)
1307
1308     fclose(param);
1309     delete [] filtmp;
1310     delete [] savemu;
1311     delete [] savesigma;
1312 }
1313 //______________________________________________________________________
1314 void AliITSsimulationSDD::Compress1D(){
1315     // 1D zero-suppression algorithm (from Gianluca A.)
1316     Int_t    dis,tol,thres,decr,diff;
1317     UChar_t *str=fStream->Stream();
1318     Int_t    counter=0;
1319     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1320
1321     Bool_t   do10to8=res->Do10to8();
1322     Int_t    last=0;
1323     Int_t    k,i,j;
1324     SetCompressParam();
1325     for (k=0; k<2; k++) {
1326         tol = Tolerance(k);
1327         dis = Disable(k);
1328         for (i=0; i<fNofMaps/2; i++) {
1329             Bool_t firstSignal=kTRUE;
1330             Int_t idx=i+k*fNofMaps/2;
1331             if( !fAnodeFire[idx] ) continue;
1332             CompressionParam(idx,decr,thres);
1333
1334             for (j=0; j<fMaxNofSamples; j++) {
1335                 Int_t signal=(Int_t)(fHitMap2->GetSignal(idx,j));
1336                 signal -= decr;  // if baseline eq.
1337                 if(do10to8) signal = Convert10to8(signal);
1338                 if (signal <= thres) {
1339                     signal=0;
1340                     diff=128;
1341                     last=0;
1342                     // write diff in the buffer for HuffT
1343                     str[counter]=(UChar_t)diff;
1344                     counter++;
1345                     continue;
1346                 } // end if signal <= thres
1347                 diff=signal-last;
1348                 if (diff > 127) diff=127;
1349                 if (diff < -128) diff=-128;
1350                 if (signal < dis) {
1351                     // tol has changed to 8 possible cases ? - one can write
1352                     // this if(TMath::Abs(diff)<tol) ... else ...
1353                     if(TMath::Abs(diff)<tol) diff=0;
1354                     // or keep it as it was before
1355                     AddDigit(idx,j,last+diff);
1356                 } else {
1357                     AddDigit(idx,j,signal);
1358                 } // end if singal < dis
1359                 diff += 128;
1360                 // write diff in the buffer used to compute Huffman tables
1361                 if (firstSignal) str[counter]=(UChar_t)signal;
1362                 else str[counter]=(UChar_t)diff;
1363                 counter++;
1364                 last=signal;
1365                 firstSignal=kFALSE;
1366             } // end for j loop time samples
1367         } // end for i loop anodes  one half of detector
1368     } //  end for k
1369
1370     // check
1371     fStream->CheckCount(counter);
1372
1373     // open file and write out the stream of diff's
1374     static Bool_t open=kTRUE;
1375     static TFile *outFile;
1376     Bool_t write = res->OutputOption();
1377     TDirectory *savedir = gDirectory;
1378
1379     if (write ) {
1380         if(open) {
1381             SetFileName("stream.root");
1382             cout<<"filename "<<fFileName<<endl;
1383             outFile=new TFile(fFileName,"recreate");
1384             cout<<"I have opened "<<fFileName<<" file "<<endl;
1385         } // end if open
1386         open = kFALSE;
1387         outFile->cd();
1388         fStream->Write();
1389     }  // endif write
1390
1391     fStream->ClearStream();
1392
1393     // back to galice.root file
1394     if(savedir) savedir->cd();
1395 }
1396 //______________________________________________________________________
1397 void AliITSsimulationSDD::StoreAllDigits(){
1398     // if non-zero-suppressed data
1399     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1400
1401     Bool_t do10to8 = res->Do10to8();
1402     Int_t i, j, digits[3];
1403
1404     for (i=0; i<fNofMaps; i++) {
1405         for (j=0; j<fMaxNofSamples; j++) {
1406             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1407             if(do10to8) signal = Convert10to8(signal);
1408             digits[0] = i;
1409             digits[1] = j;
1410             digits[2] = signal;
1411             fITS->AddRealDigit(1,digits);
1412         } // end for j
1413     } // end for i
1414 }
1415 //______________________________________________________________________
1416 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
1417     // Creates histograms of maps for debugging
1418     Int_t i;
1419
1420     fHis=new TObjArray(fNofMaps);
1421     for (i=0;i<fNofMaps;i++) {
1422         TString sddName("sdd_");
1423         Char_t candNum[4];
1424         sprintf(candNum,"%d",i+1);
1425         sddName.Append(candNum);
1426         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
1427                              0.,(Float_t) scale*fMaxNofSamples), i);
1428     } // end for i
1429 }
1430 //______________________________________________________________________
1431 void AliITSsimulationSDD::FillHistograms(){
1432     // fill 1D histograms from map
1433
1434     if (!fHis) return;
1435
1436     for( Int_t i=0; i<fNofMaps; i++) {
1437         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
1438         Int_t nsamples = hist->GetNbinsX();
1439         for( Int_t j=0; j<nsamples; j++) {
1440             Double_t signal=fHitMap2->GetSignal(i,j);
1441             hist->Fill((Float_t)j,signal);
1442         } // end for j
1443     } // end for i
1444 }
1445 //______________________________________________________________________
1446 void AliITSsimulationSDD::ResetHistograms(){
1447     // Reset histograms for this detector
1448     Int_t i;
1449
1450     for (i=0;i<fNofMaps;i++ ) {
1451         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
1452     } // end for i
1453 }
1454 //______________________________________________________________________
1455 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
1456     // Fills a histogram from a give anode.
1457
1458     if (!fHis) return 0;
1459
1460     if(wing <=0 || wing > 2) {
1461         Warning("GetAnode","Wrong wing number: %d",wing);
1462         return NULL;
1463     } // end if wing <=0 || wing >2
1464     if(anode <=0 || anode > fNofMaps/2) {
1465         Warning("GetAnode","Wrong anode number: %d",anode);
1466         return NULL;
1467     } // end if ampde <=0 || andoe > fNofMaps/2
1468
1469     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
1470     return (TH1F*)(fHis->At(index));
1471 }
1472 //______________________________________________________________________
1473 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
1474     // Writes the histograms to a file
1475
1476     if (!fHis) return;
1477
1478     hfile->cd();
1479     Int_t i;
1480     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
1481     return;
1482 }
1483 //______________________________________________________________________
1484 Float_t AliITSsimulationSDD::GetNoise() {
1485     // Returns the noise value
1486     //Bool_t do10to8=GetResp()->Do10to8();
1487     //noise will always be in the liniar part of the signal
1488     Int_t decr;
1489     Int_t threshold = fT1[0];
1490     char opt1[20], opt2[20];
1491     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1492     SetCompressParam();
1493     res->GetParamOptions(opt1,opt2);
1494     fParam=opt2;
1495     Double_t noise,baseline;
1496     //GetBaseline(fModule);
1497
1498     TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
1499     if(c2) delete c2->GetPrimitive("noisehist");
1500     if(c2) delete c2->GetPrimitive("anode");
1501     else     c2=new TCanvas("c2");
1502     c2->cd();
1503     c2->SetFillColor(0);
1504
1505     TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold);
1506     TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0.,
1507                            (float)fMaxNofSamples);
1508     Int_t i,k;
1509     for (i=0;i<fNofMaps;i++) {
1510         CompressionParam(i,decr,threshold);
1511         baseline = res->GetBaseline(i);
1512         noise = res->GetNoise(i);
1513         anode->Reset();
1514         for (k=0;k<fMaxNofSamples;k++) {
1515             Float_t signal=(Float_t)fHitMap2->GetSignal(i,k);
1516             //if (signal <= (float)threshold) noisehist->Fill(signal-baseline);
1517             if (signal <= (float)(threshold+decr)) noisehist->Fill(signal);
1518             anode->Fill((float)k,signal);
1519         } // end for k
1520         anode->Draw();
1521         c2->Update();
1522     } // end for i
1523     TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold);
1524     noisehist->Fit("gnoise","RQ");
1525     noisehist->Draw();
1526     c2->Update();
1527     Float_t mnoise = gnoise->GetParameter(1);
1528     cout << "mnoise : " << mnoise << endl;
1529     Float_t rnoise = gnoise->GetParameter(2);
1530     cout << "rnoise : " << rnoise << endl;
1531     delete noisehist;
1532     return rnoise;
1533 }
1534 //______________________________________________________________________
1535 void AliITSsimulationSDD::WriteSDigits(){
1536     // Fills the Summable digits Tree
1537     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1538
1539     for( Int_t i=0; i<fNofMaps; i++ ) {
1540         if( !fAnodeFire[i] ) continue;
1541         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1542             Double_t sig = fHitMap2->GetSignal( i, j );
1543             if( sig > 0.2 ) {
1544                 Int_t jdx = j*fScaleSize;
1545                 Int_t index = fpList->GetHitIndex( i, j );
1546                 AliITSpListItem pItemTmp2( fModule, index, 0. );
1547                 // put the fScaleSize analog digits in only one
1548                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1549                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1550                     if( pItemTmp == 0 ) continue;
1551                     pItemTmp2.Add( pItemTmp );
1552                 }
1553                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1554                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1555                 aliITS->AddSumDigit( pItemTmp2 );
1556             } // end if (sig > 0.2)
1557         }
1558     }
1559     return;
1560 }
1561 //______________________________________________________________________
1562 void AliITSsimulationSDD::PrintStatus() const {
1563     // Print SDD simulation Parameters
1564
1565     cout << "**************************************************" << endl;
1566     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1567     cout << "**************************************************" << endl;
1568     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1569     cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl;
1570     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1571     cout << "Number pf Anodes used: " << fNofMaps << endl;
1572     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1573     cout << "Scale size factor: " << fScaleSize << endl;
1574     cout << "**************************************************" << endl;
1575 }