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