ITS/AliITSsimulationSDD.cxx

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