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)) ){
 754       printf("Chip bad mod %d  chip %d anode %d\n",mod,res->GetChip(i),i);
 755       gain=0.;
 756     }
 757     for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 758       fInZR[k]  = fHitMap2->GetSignal(i,k);
 759       if(bAddGain) fInZR[k]*=gain;
 760       if( bAddNoise ) {
 761         contrib   = (baseline + noise*gRandom->Gaus());
 762         fInZR[k] += contrib;
 763       }
 764       fInZI[k]  = 0.;
 765     } // end for k
 766     if(!fDoFFT) {
 767       for(k=0; k<fMaxNofSamples; k++) {
 768         Double_t newcont = 0.;
 769         Double_t maxcont = 0.;
 770         for(kk=0;kk<fScaleSize;kk++) {
 771           newcont = fInZR[fScaleSize*k+kk];
 772           if(newcont > maxcont) maxcont = newcont;
 773         } // end for kk
 774         newcont = maxcont;
 775         if (newcont >= maxadc) newcont = maxadc -1;
 776         if(newcont >= baseline){
 777           Warning("","newcont=%d>=baseline=%d",newcont,baseline);
 778         } // end if
 779           // back to analog: ?
 780         fHitMap2->SetHit(i,k,newcont);
 781       }  // end for k
 782     }else{
 783       FastFourierTransform(fElectronics,&fInZR[0],&fInZI[0],1);
 784       for(k=0; k<fScaleSize*fMaxNofSamples; k++) {
 785         Double_t rw = fElectronics->GetTraFunReal(k);
 786         Double_t iw = fElectronics->GetTraFunImag(k);
 787         fOutZR[k]   = fInZR[k]*rw - fInZI[k]*iw;
 788         fOutZI[k]   = fInZR[k]*iw + fInZI[k]*rw;
 789       } // end for k
 790       FastFourierTransform(fElectronics,&fOutZR[0],&fOutZI[0],-1);
 791       for(k=0; k<fMaxNofSamples; k++) {
 792         Double_t newcont1 = 0.;
 793         Double_t maxcont1 = 0.;
 794         for(kk=0;kk<fScaleSize;kk++) {
 795           newcont1 = fOutZR[fScaleSize*k+kk];
 796           if(newcont1 > maxcont1) maxcont1 = newcont1;
 797         } // end for kk
 798         newcont1 = maxcont1;
 799         if (newcont1 >= maxadc) newcont1 = maxadc -1;
 800         fHitMap2->SetHit(i,k,newcont1);
 801       } // end for k
 802     }
 803   } // end for i loop over anodes
 804   return;
 805 }
 806
 807 //______________________________________________________________________
 808 void AliITSsimulationSDD::ApplyCrosstalk(Int_t mod) {
 809     // function add the crosstalk effect to signal
 810     // temporal function, should be checked...!!!
 811     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
 812
 813     Int_t fNofMaps = seg->Npz();
 814     Int_t fMaxNofSamples = seg->Npx();
 815
 816     // create and inizialice crosstalk map
 817     Float_t* ctk = new Float_t[fNofMaps*fMaxNofSamples+1];
 818     if( ctk == NULL ) {
 819         Error( "ApplyCrosstalk", "no memory for temporal map: exit \n" );
 820         return;
 821     }
 822     memset( ctk, 0, sizeof(Float_t)*(fNofMaps*fMaxNofSamples+1) );
 823     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(mod);
 824     for( Int_t z=0; z<fNofMaps; z++ ) {
 825       Double_t baseline = calibr->GetBaseline(z);
 826         Bool_t on = kFALSE;
 827         Int_t tstart = 0;
 828         Int_t tstop = 0;
 829         Int_t nTsteps = 0;
 830
 831         for( Int_t l=0; l<fMaxNofSamples; l++ ) {
 832             Float_t fadc = (Float_t)fHitMap2->GetSignal( z, l );
 833             if( fadc > baseline ) {
 834                 if( on == kFALSE && l<fMaxNofSamples-4 ) {
 835                     Float_t fadc1 = (Float_t)fHitMap2->GetSignal( z, l+1 );
 836                     if( fadc1 < fadc ) continue;
 837                     on = kTRUE;
 838                     nTsteps = 0;
 839                     tstart = l;
 840                 }
 841                 nTsteps++;
 842             }
 843             else { // end fadc > baseline
 844                 if( on == kTRUE ) {
 845                     if( nTsteps > 2 ) {
 846                         tstop = l;
 847                         // make smooth derivative
 848                         Float_t* dev = new Float_t[fMaxNofSamples+1];
 849                         memset( dev, 0, sizeof(Float_t)*(fMaxNofSamples+1) );
 850                         if( ctk == NULL ) {
 851                             Error( "ApplyCrosstalk",
 852                                    "no memory for temporal array: exit \n" );
 853                             return;
 854                         }
 855                         for( Int_t i=tstart; i<tstop; i++ ) {
 856                             if( i > 2 && i < fMaxNofSamples-2 )
 857                                 dev[i] = -0.2*fHitMap2->GetSignal( z,i-2 )
 858                                     -0.1*fHitMap2->GetSignal( z,i-1 )
 859                                     +0.1*fHitMap2->GetSignal( z,i+1 )
 860                                     +0.2*fHitMap2->GetSignal( z,i+2 );
 861                         }
 862
 863                         // add crosstalk contribution to neibourg anodes
 864                         for( Int_t i=tstart; i<tstop; i++ ) {
 865                             Int_t anode = z - 1;
 866                             Int_t i1 = (Int_t)((i-tstart)*.61+tstart+0.5); //
 867                             Float_t ctktmp =  -dev[i1] * 0.25;
 868                             if( anode > 0 ) {
 869                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 870                             }
 871                             anode = z + 1;
 872                             if( anode < fNofMaps ) {
 873                                 ctk[anode*fMaxNofSamples+i] += ctktmp;
 874                             }
 875                         }
 876                         delete [] dev;
 877
 878                     } // if( nTsteps > 2 )
 879                     on = kFALSE;
 880                 }  // if( on == kTRUE )
 881             }  // else
 882         }
 883     }
 884
 885     for( Int_t a=0; a<fNofMaps; a++ )
 886         for( Int_t t=0; t<fMaxNofSamples; t++ ) {
 887             Float_t signal = fHitMap2->GetSignal(a,t)+ctk[a*fMaxNofSamples+t];
 888             fHitMap2->SetHit( a, t, signal );
 889         }
 890
 891     delete [] ctk;
 892 }
 893
 894 //______________________________________________________________________
 895 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl,
 896                                            Int_t &th) const{
 897     // Returns the compression alogirthm parameters
 898   db=fD[i];
 899   tl=fT1[i];
 900   th=fT2[i];
 901 }
 902 //______________________________________________________________________
 903 void AliITSsimulationSDD::CompressionParam(Int_t i,Int_t &db,Int_t &tl) const{
 904     // returns the compression alogirthm parameters
 905
 906     db=fD[i];
 907     tl=fT1[i];
 908
 909 }
 910 //______________________________________________________________________
 911 void AliITSsimulationSDD::SetCompressParam(){
 912     // Sets the compression alogirthm parameters
 913     AliITSCalibrationSDD* calibr = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 914     for(Int_t ian = 0; ian<fNofMaps;ian++){
 915       fD[ian] = (Int_t)(calibr->GetBaseline(ian));
 916       fT1[ian] = (Int_t)(2.*calibr->GetNoiseAfterElectronics(ian)+0.5);
 917       fT2[ian] = 0;   // used by 2D clustering - not defined yet
 918       fTol[ian] = 0; // used by 2D clustering - not defined yet
 919     }
 920 }
 921 //______________________________________________________________________
 922 Int_t AliITSsimulationSDD::Convert10to8(Int_t signal) const {
 923     // To the 10 to 8 bit lossive compression.
 924     // code from Davide C. and Albert W.
 925
 926     if (signal < 128)  return signal;
 927     if (signal < 256)  return (128+((signal-128)>>1));
 928     if (signal < 512)  return (192+((signal-256)>>3));
 929     if (signal < 1024) return (224+((signal-512)>>4));
 930     return 0;
 931 }
 932 /*
 933 //______________________________________________________________________
 934 AliITSMap*   AliITSsimulationSDD::HitMap(Int_t i){
 935     //Return the correct map.
 936
 937     return ((i==0)? fHitMap1 : fHitMap2);
 938 }
 939 */
 940 //______________________________________________________________________
 941 void AliITSsimulationSDD::ZeroSuppression(const char *option) {
 942     // perform the zero suppresion
 943     if (strstr(option,"2D")) {
 944         //Init2D();              // activate if param change module by module
 945         Compress2D();
 946     } else if (strstr(option,"1D")) {
 947         //Init1D();              // activate if param change module by module
 948         Compress1D();
 949     } else StoreAllDigits();
 950 }
 951 //______________________________________________________________________
 952 void AliITSsimulationSDD::Init2D(){
 953     // read in and prepare arrays: fD, fT1, fT2
 954     //                         savemu[nanodes], savesigma[nanodes]
 955     // read baseline and noise from file - either a .root file and in this
 956     // case data should be organised in a tree with one entry for each
 957     // module => reading should be done accordingly
 958     // or a classic file and do smth. like this ( code from Davide C. and
 959     // Albert W.) :
 960     // Read 2D zero-suppression parameters for SDD
 961
 962     if (!strstr(fParam.Data(),"file")) return;
 963
 964     Int_t na,pos,tempTh;
 965     Float_t mu,sigma;
 966     Float_t *savemu    = new Float_t [fNofMaps];
 967     Float_t *savesigma = new Float_t [fNofMaps];
 968     char input[100],basel[100],par[100];
 969     char *filtmp;
 970     Double_t tmp1,tmp2;
 971     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
 972
 973     res->Thresholds(tmp1,tmp2);
 974     Int_t minval = static_cast<Int_t>(tmp1);
 975
 976     res->Filenames(input,basel,par);
 977     fFileName = par;
 978     //
 979     filtmp = gSystem->ExpandPathName(fFileName.Data());
 980     FILE *param = fopen(filtmp,"r");
 981     na = 0;
 982
 983     if(param) {
 984         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
 985             if (pos != na+1) {
 986                 Error("Init2D","Anode number not in increasing order!",filtmp);
 987                 exit(1);
 988             } // end if pos != na+1
 989             savemu[na] = mu;
 990             savesigma[na] = sigma;
 991             if ((2.*sigma) < mu) {
 992                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
 993                 mu = 2.0 * sigma;
 994             } else fD[na] = 0;
 995             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
 996             if (tempTh < 0) tempTh=0;
 997             fT1[na] = tempTh;
 998             tempTh = (Int_t)floor(mu+3.0*sigma+0.5) - minval;
 999             if (tempTh < 0) tempTh=0;
1000             fT2[na] = tempTh;
1001             na++;
1002         } // end while
1003     } else {
1004         Error("Init2D","THE FILE %s DOES NOT EXIST !",filtmp);
1005         exit(1);
1006     } // end if(param)
1007
1008     fclose(param);
1009     delete [] filtmp;
1010     delete [] savemu;
1011     delete [] savesigma;
1012 }
1013 //______________________________________________________________________
1014 void AliITSsimulationSDD::Compress2D(){
1015     // simple ITS cluster finder -- online zero-suppression conditions
1016
1017     Int_t db,tl,th;
1018     Double_t tmp1,tmp2;
1019     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1020
1021     res->Thresholds(tmp1,tmp2);
1022     Int_t minval   = static_cast<Int_t>(tmp1);
1023     Bool_t write   = res->OutputOption();
1024     Bool_t do10to8 = res->Do10to8();
1025     Int_t nz, nl, nh, low, i, j;
1026     SetCompressParam();
1027     for (i=0; i<fNofMaps; i++) {
1028         CompressionParam(i,db,tl,th);
1029         nz  = 0;
1030         nl  = 0;
1031         nh  = 0;
1032         low = 0;
1033         for (j=0; j<fMaxNofSamples; j++) {
1034             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1035             signal -= db; // if baseline eq. is done here
1036             if (signal <= 0) {nz++; continue;}
1037             if ((signal - tl) < minval) low++;
1038             if ((signal - th) >= minval) {
1039                 nh++;
1040                 Bool_t cond=kTRUE;
1041                 FindCluster(i,j,signal,minval,cond);
1042                 if(cond && j &&
1043                    ((TMath::Abs(fHitMap2->GetSignal(i,j-1))-th)>=minval)){
1044                     if(do10to8) signal = Convert10to8(signal);
1045                     AddDigit(i,j,signal);
1046                 } // end if cond&&j&&()
1047             } else if ((signal - tl) >= minval) nl++;
1048         } // end for j loop time samples
1049         if (write) TreeB()->Fill(nz,nl,nh,low,i+1);
1050     } //end for i loop anodes
1051
1052     char hname[30];
1053     if (write) {
1054         sprintf(hname,"TNtuple%d_%d",fModule,fEvent);
1055         TreeB()->Write(hname);
1056         // reset tree
1057         TreeB()->Reset();
1058     } // end if write
1059 }
1060 //______________________________________________________________________
1061 void  AliITSsimulationSDD::FindCluster(Int_t i,Int_t j,Int_t signal,
1062                                        Int_t minval,Bool_t &cond){
1063     // Find clusters according to the online 2D zero-suppression algorithm
1064     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1065     AliITSsegmentationSDD* seg = (AliITSsegmentationSDD*)GetSegmentationModel(1);
1066
1067     Bool_t do10to8 = res->Do10to8();
1068     Bool_t high    = kFALSE;
1069
1070     fHitMap2->FlagHit(i,j);
1071     //
1072     //  check the online zero-suppression conditions
1073     //
1074     const Int_t kMaxNeighbours = 4;
1075     Int_t nn;
1076     Int_t dbx,tlx,thx;
1077     Int_t xList[kMaxNeighbours], yList[kMaxNeighbours];
1078     seg->Neighbours(i,j,&nn,xList,yList);
1079     Int_t in,ix,iy,qns;
1080     for (in=0; in<nn; in++) {
1081         ix=xList[in];
1082         iy=yList[in];
1083         if (fHitMap2->TestHit(ix,iy)==kUnused) {
1084             CompressionParam(ix,dbx,tlx,thx);
1085             Int_t qn = (Int_t)(fHitMap2->GetSignal(ix,iy));
1086             qn -= dbx; // if baseline eq. is done here
1087             if ((qn-tlx) < minval) {
1088                 fHitMap2->FlagHit(ix,iy);
1089                 continue;
1090             } else {
1091                 if ((qn - thx) >= minval) high=kTRUE;
1092                 if (cond) {
1093                     if(do10to8) signal = Convert10to8(signal);
1094                     AddDigit(i,j,signal);
1095                 } // end if cond
1096                 if(do10to8) qns = Convert10to8(qn);
1097                 else qns=qn;
1098                 if (!high) AddDigit(ix,iy,qns);
1099                 cond=kFALSE;
1100                 if(!high) fHitMap2->FlagHit(ix,iy);
1101             } // end if qn-tlx < minval
1102         } // end if  TestHit
1103     } // end for in loop over neighbours
1104 }
1105 //______________________________________________________________________
1106 void AliITSsimulationSDD::Init1D(){
1107     // this is just a copy-paste of input taken from 2D algo
1108     // Torino people should give input
1109     // Read 1D zero-suppression parameters for SDD
1110
1111     if (!strstr(fParam.Data(),"file")) return;
1112
1113     Int_t na,pos,tempTh;
1114     Float_t mu,sigma;
1115     Float_t *savemu    = new Float_t [fNofMaps];
1116     Float_t *savesigma = new Float_t [fNofMaps];
1117     char input[100],basel[100],par[100];
1118     char *filtmp;
1119     Double_t tmp1,tmp2;
1120     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1121
1122     res->Thresholds(tmp1,tmp2);
1123     Int_t minval = static_cast<Int_t>(tmp1);
1124
1125     res->Filenames(input,basel,par);
1126     fFileName=par;
1127
1128     //  set first the disable and tol param
1129     SetCompressParam();
1130     //
1131     filtmp = gSystem->ExpandPathName(fFileName.Data());
1132     FILE *param = fopen(filtmp,"r");
1133     na = 0;
1134
1135     if (param) {
1136         fscanf(param,"%d %d %d %d ", &fT2[0], &fT2[1], &fTol[0], &fTol[1]);
1137         while(fscanf(param,"%d %f %f",&pos, &mu, &sigma) != EOF) {
1138             if (pos != na+1) {
1139                 Error("Init1D","Anode number not in increasing order!",filtmp);
1140                 exit(1);
1141             } // end if pos != na+1
1142             savemu[na]=mu;
1143             savesigma[na]=sigma;
1144             if ((2.*sigma) < mu) {
1145                 fD[na] = (Int_t)floor(mu - 2.0*sigma + 0.5);
1146                 mu = 2.0 * sigma;
1147             } else fD[na] = 0;
1148             tempTh = (Int_t)floor(mu+2.25*sigma+0.5) - minval;
1149             if (tempTh < 0) tempTh=0;
1150             fT1[na] = tempTh;
1151             na++;
1152         } // end while
1153     } else {
1154         Error("Init1D","THE FILE %s DOES NOT EXIST !",filtmp);
1155         exit(1);
1156     } // end if(param)
1157
1158     fclose(param);
1159     delete [] filtmp;
1160     delete [] savemu;
1161     delete [] savesigma;
1162 }
1163 //______________________________________________________________________
1164 void AliITSsimulationSDD::Compress1D(){
1165     // 1D zero-suppression algorithm (from Gianluca A.)
1166     Int_t    dis,tol,thres,decr,diff;
1167     UChar_t *str=fStream->Stream();
1168     Int_t    counter=0;
1169     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1170
1171     Bool_t   do10to8=res->Do10to8();
1172     Int_t    last=0;
1173     Int_t    k,i,j;
1174     SetCompressParam();
1175     for (k=0; k<2; k++) {
1176         tol = Tolerance(k);
1177         dis = Disable(k);
1178         for (i=0; i<fNofMaps/2; i++) {
1179             Bool_t firstSignal=kTRUE;
1180             Int_t idx=i+k*fNofMaps/2;
1181             if( !fAnodeFire[idx] ) continue;
1182             CompressionParam(idx,decr,thres);
1183
1184             for (j=0; j<fMaxNofSamples; j++) {
1185                 Int_t signal=(Int_t)(fHitMap2->GetSignal(idx,j));
1186                 signal -= decr;  // if baseline eq.
1187                 if(do10to8) signal = Convert10to8(signal);
1188                 if (signal <= thres) {
1189                     signal=0;
1190                     diff=128;
1191                     last=0;
1192                     // write diff in the buffer for HuffT
1193                     str[counter]=(UChar_t)diff;
1194                     counter++;
1195                     continue;
1196                 } // end if signal <= thres
1197                 diff=signal-last;
1198                 if (diff > 127) diff=127;
1199                 if (diff < -128) diff=-128;
1200                 if (signal < dis) {
1201                     // tol has changed to 8 possible cases ? - one can write
1202                     // this if(TMath::Abs(diff)<tol) ... else ...
1203                     if(TMath::Abs(diff)<tol) diff=0;
1204                     // or keep it as it was before
1205                     AddDigit(idx,j,last+diff);
1206                 } else {
1207                     AddDigit(idx,j,signal);
1208                 } // end if singal < dis
1209                 diff += 128;
1210                 // write diff in the buffer used to compute Huffman tables
1211                 if (firstSignal) str[counter]=(UChar_t)signal;
1212                 else str[counter]=(UChar_t)diff;
1213                 counter++;
1214                 last=signal;
1215                 firstSignal=kFALSE;
1216             } // end for j loop time samples
1217         } // end for i loop anodes  one half of detector
1218     } //  end for k
1219
1220     // check
1221     fStream->CheckCount(counter);
1222
1223     // open file and write out the stream of diff's
1224     static Bool_t open=kTRUE;
1225     static TFile *outFile;
1226     Bool_t write = res->OutputOption();
1227     TDirectory *savedir = gDirectory;
1228
1229     if (write ) {
1230         if(open) {
1231             SetFileName("stream.root");
1232             cout<<"filename "<<fFileName<<endl;
1233             outFile=new TFile(fFileName,"recreate");
1234             cout<<"I have opened "<<fFileName<<" file "<<endl;
1235         } // end if open
1236         open = kFALSE;
1237         outFile->cd();
1238         fStream->Write();
1239     }  // endif write
1240
1241     fStream->ClearStream();
1242
1243     // back to galice.root file
1244     if(savedir) savedir->cd();
1245 }
1246 //______________________________________________________________________
1247 void AliITSsimulationSDD::StoreAllDigits(){
1248     // if non-zero-suppressed data
1249     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1250
1251     Bool_t do10to8 = res->Do10to8();
1252     Int_t i, j, digits[3];
1253
1254     for (i=0; i<fNofMaps; i++) {
1255         for (j=0; j<fMaxNofSamples; j++) {
1256             Int_t signal=(Int_t)(fHitMap2->GetSignal(i,j));
1257             if(do10to8) signal = Convert10to8(signal);
1258             digits[0] = i;
1259             digits[1] = j;
1260             digits[2] = signal;
1261             fITS->AddRealDigit(1,digits);
1262         } // end for j
1263     } // end for i
1264 }
1265 //______________________________________________________________________
1266 void AliITSsimulationSDD::CreateHistograms(Int_t scale){
1267     // Creates histograms of maps for debugging
1268     Int_t i;
1269
1270     fHis=new TObjArray(fNofMaps);
1271     for (i=0;i<fNofMaps;i++) {
1272         TString sddName("sdd_");
1273         Char_t candNum[4];
1274         sprintf(candNum,"%d",i+1);
1275         sddName.Append(candNum);
1276         fHis->AddAt(new TH1F(sddName.Data(),"SDD maps",scale*fMaxNofSamples,
1277                              0.,(Float_t) scale*fMaxNofSamples), i);
1278     } // end for i
1279 }
1280 //______________________________________________________________________
1281 void AliITSsimulationSDD::FillHistograms(){
1282     // fill 1D histograms from map
1283
1284     if (!fHis) return;
1285
1286     for( Int_t i=0; i<fNofMaps; i++) {
1287         TH1F *hist =(TH1F *)fHis->UncheckedAt(i);
1288         Int_t nsamples = hist->GetNbinsX();
1289         for( Int_t j=0; j<nsamples; j++) {
1290             Double_t signal=fHitMap2->GetSignal(i,j);
1291             hist->Fill((Float_t)j,signal);
1292         } // end for j
1293     } // end for i
1294 }
1295 //______________________________________________________________________
1296 void AliITSsimulationSDD::ResetHistograms(){
1297     // Reset histograms for this detector
1298     Int_t i;
1299
1300     for (i=0;i<fNofMaps;i++ ) {
1301         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
1302     } // end for i
1303 }
1304 //______________________________________________________________________
1305 TH1F *AliITSsimulationSDD::GetAnode(Int_t wing, Int_t anode) {
1306     // Fills a histogram from a give anode.
1307
1308     if (!fHis) return 0;
1309
1310     if(wing <=0 || wing > 2) {
1311         Warning("GetAnode","Wrong wing number: %d",wing);
1312         return NULL;
1313     } // end if wing <=0 || wing >2
1314     if(anode <=0 || anode > fNofMaps/2) {
1315         Warning("GetAnode","Wrong anode number: %d",anode);
1316         return NULL;
1317     } // end if ampde <=0 || andoe > fNofMaps/2
1318
1319     Int_t index = (wing-1)*fNofMaps/2 + anode-1;
1320     return (TH1F*)(fHis->At(index));
1321 }
1322 //______________________________________________________________________
1323 void AliITSsimulationSDD::WriteToFile(TFile *hfile) {
1324     // Writes the histograms to a file
1325
1326     if (!fHis) return;
1327
1328     hfile->cd();
1329     Int_t i;
1330     for(i=0; i<fNofMaps; i++)  fHis->At(i)->Write(); //fAdcs[i]->Write();
1331     return;
1332 }
1333 //______________________________________________________________________
1334 Float_t AliITSsimulationSDD::GetNoise() {
1335     // Returns the noise value
1336     //Bool_t do10to8=GetResp()->Do10to8();
1337     //noise will always be in the liniar part of the signal
1338     Int_t decr;
1339     Int_t threshold = fT1[0];
1340     char opt1[20], opt2[20];
1341     AliITSCalibrationSDD* res = (AliITSCalibrationSDD*)GetCalibrationModel(fModule);
1342     SetCompressParam();
1343     res->GetParamOptions(opt1,opt2);
1344     fParam=opt2;
1345     Double_t noise,baseline;
1346     //GetBaseline(fModule);
1347
1348     TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("c2");
1349     if(c2) delete c2->GetPrimitive("noisehist");
1350     if(c2) delete c2->GetPrimitive("anode");
1351     else     c2=new TCanvas("c2");
1352     c2->cd();
1353     c2->SetFillColor(0);
1354
1355     TH1F *noisehist = new TH1F("noisehist","noise",100,0.,(float)2*threshold);
1356     TH1F *anode = new TH1F("anode","Anode Projection",fMaxNofSamples,0.,
1357                            (float)fMaxNofSamples);
1358     Int_t i,k;
1359     for (i=0;i<fNofMaps;i++) {
1360         CompressionParam(i,decr,threshold);
1361         baseline = res->GetBaseline(i);
1362         noise = res->GetNoise(i);
1363         anode->Reset();
1364         for (k=0;k<fMaxNofSamples;k++) {
1365             Float_t signal=(Float_t)fHitMap2->GetSignal(i,k);
1366             //if (signal <= (float)threshold) noisehist->Fill(signal-baseline);
1367             if (signal <= (float)(threshold+decr)) noisehist->Fill(signal);
1368             anode->Fill((float)k,signal);
1369         } // end for k
1370         anode->Draw();
1371         c2->Update();
1372     } // end for i
1373     TF1 *gnoise = new TF1("gnoise","gaus",0.,threshold);
1374     noisehist->Fit("gnoise","RQ");
1375     noisehist->Draw();
1376     c2->Update();
1377     Float_t mnoise = gnoise->GetParameter(1);
1378     cout << "mnoise : " << mnoise << endl;
1379     Float_t rnoise = gnoise->GetParameter(2);
1380     cout << "rnoise : " << rnoise << endl;
1381     delete noisehist;
1382     return rnoise;
1383 }
1384 //______________________________________________________________________
1385 void AliITSsimulationSDD::WriteSDigits(){
1386     // Fills the Summable digits Tree
1387     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
1388
1389     for( Int_t i=0; i<fNofMaps; i++ ) {
1390         if( !fAnodeFire[i] ) continue;
1391         for( Int_t j=0; j<fMaxNofSamples; j++ ) {
1392             Double_t sig = fHitMap2->GetSignal( i, j );
1393             if( sig > 0.2 ) {
1394                 Int_t jdx = j*fScaleSize;
1395                 Int_t index = fpList->GetHitIndex( i, j );
1396                 AliITSpListItem pItemTmp2( fModule, index, 0. );
1397                 // put the fScaleSize analog digits in only one
1398                 for( Int_t ik=0; ik<fScaleSize; ik++ ) {
1399                     AliITSpListItem *pItemTmp = fpList->GetpListItem(i,jdx+ik);
1400                     if( pItemTmp == 0 ) continue;
1401                     pItemTmp2.Add( pItemTmp );
1402                 }
1403                 pItemTmp2.AddSignalAfterElect( fModule, index, sig );
1404                 pItemTmp2.AddNoise(fModule,index,fHitNoiMap2->GetSignal(i,j));
1405                 aliITS->AddSumDigit( pItemTmp2 );
1406             } // end if (sig > 0.2)
1407         }
1408     }
1409     return;
1410 }
1411 //______________________________________________________________________
1412 void AliITSsimulationSDD::PrintStatus() const {
1413     // Print SDD simulation Parameters
1414
1415     cout << "**************************************************" << endl;
1416     cout << "   Silicon Drift Detector Simulation Parameters   " << endl;
1417     cout << "**************************************************" << endl;
1418     cout << "Flag for Perpendicular tracks: " << (Int_t) fFlag << endl;
1419     cout << "Flag for noise checking: " << (Int_t) fCheckNoise << endl;
1420     cout << "Flag to switch off electronics: " << (Int_t) fDoFFT << endl;
1421     cout << "Number pf Anodes used: " << fNofMaps << endl;
1422     cout << "Number of Time Samples: " << fMaxNofSamples << endl;
1423     cout << "Scale size factor: " << fScaleSize << endl;
1424     cout << "**************************************************" << endl;
1425 }