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