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